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Abstract 

Numerical solution of infinite-domain boundary- 
value problems requires some special techniques that 
would make the problem available for treatment on 
the computer. Indeed, the problem must be dis- 
cretized in a way that the computer operates with 
only finite amount of information. Therefore, the 
original infinite-domain formulation must be altered 
and/or augmented so that on one hand the solution 
is not changed (or changed slightly) and on the other 
hand the finite discrete formulation becomes avail- 
able. 

One widely used approach to constructing such 
discretizations consists of truncating the unbounded 
original domain and then setting the artificial 
boundary conditions (ABC’s) at the newly formed 
external boundary. The role of the ABC’s is to close 
the truncated problem and at the same time to en- 
sure that the solution found inside the finite compu- 
tational domain would be maximally close to (in the 
ideal case, exactly the same as) the corresponding 
fragment of the original infinite-domain solution. 

Let us emphasize that the proper treatment of ar- 
tificial boundaries may have a profound impact on 
the overall quality and performance of numerical al- 
gorithms. The latter statement is corroborated by 
the numerous computational experiments and espe- 
cially concerns the area of CFD, in which external 
problems present a wide class of practically impor- 
tant formulations. 

In this paper, we review some work that has been 
done over the recent years on constructing highly ac- 
curate nonlocal ABC’s for calculation of compress- 
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ible external flows. The approach is based on im- 
plementation of the generalized potentials and pseu- 
dodifferential boundary projection operators analo- 
gous to those proposed first by Calderon. The dif- 
ference potentials method (DPM) by Ryaben’kii is 
used for the effective computation of the generalized 
potentials and projections. The resulting ABC’s 
clearly outperform the existing methods from the 
standpoints of accuracy and robustness, in many 
cases noticeably speed up the multigrid convergence, 
and at the same time are quite comparable to other 
methods from the standpoints of geometric univer- 
sality and easiness of implementation. 

1 Basic Ideas 
1.1 Preliminaries 

Let us consider two domains: the finite computa- 
tional domain Di n and its infinite complement to 
the entire plane (entire space) D ex ; the domains 
are separated by the artificial boundary T, which 
is supposed to be a simple closed curve (surface). 
We originally formulate our problem on Di n U D ex \ 
specifically, we are looking for a vector-function 
u = u(r, t), r represents the space coordinates and 
t is the time, that satisfies some (generally speak- 
ing, nonlinear) system of partial differential equa- 
tions (PDE’s) 

$ ( U) Vu >---) = S( r > 0> ( la ) 

r £ Din U D ex , 

as well as some boundary conditions at infinity 

<p(u) — » 0, as r = |r| — *■ oo. (lb) 

We will always assume that the problem (1) is 
uniquely solvable and well-posed. Later on, we 
will think of u as of the hydrodynamic variables, 
e.g., velocity components, density, and pressure. In 
this case, system (la) is one of the relevant sys- 
tems of PDE’s used in fluid dynamics, e.g., the 
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Navier-Stokes equations, and boundary conditions 
(lb) may reflect, e.g., the fact that as the coor- 
dinate approaches infinity the local flow parame- 
ters approach the free-stream ones. We will also 
assume that in the case of an external flow prob- 
lem the computational domain Di n entirely contains 
the immersed body(ies), and equations (la) include 
the proper boundary conditions at the solid sur- 
face^). In this case, we may think (for simplicity) 
that Di n U D ex = R 11 , n = 2 or n = 3 . 

The next step is to assume that in the domain D ex 
system (la) is either linear and homogeneous from 
the very beginning or admits replacement by some 
linear homogeneous system of PDF’s under certain 
conditions. Analogously, we will assume that bound- 
ary conditions (lb) can also be thought of as some 
linear homogeneous relations. Then, instead of (1) 
we consider 

$ ( u > Vu > "•) = S( r > *)> (2a) 


Before proceeding further, let us define here an im- 
portant concept of exact ABC’s, which we will refer 
to henceforth. Namely, if the ABC’s (3b) are con- 
structed so that the solution to (3) and (2) coincide 
on Di n , then we will call such boundary conditions 
the exact ABC’s. Note, we always assume that the 
solutions of (1) and (2) are close enough and always 
treat the exactness of ABC’s within the accuracy of 
replacement of (1) by (2) (e.g., within the accuracy 
of far-field linearization). Clearly, we can reformu- 
late the definition of exact ABC’s as follows: one 
should be able to uniquely complement the solution 
of (3) from Di n to D ex so that the resulting function 
on D^ U Dex solves (2). 

An extensive work has been done by many authors 
over the recent years towards constructing the exact 
ABC’s for different problems, as well as towards de- 
veloping the effective approximate boundary condi- 
tions. A survey of this work can be found, e.g., in 
the recent review paper by Tsynkov 1 , as well as in 
the comprehensive reviews by Givoli . 


Lu = 0, r eD ex , (2b) 

lu — > 0, as r — r oo. (2c) 

In (2b) and (2c), L and 1 are some linear operators, 
and u = u(r, t ) may either coincide with u or be 
some function of u. For example, system (2b) can 
be obtained on the basis of (la) using linearization in 
D ex . (Note, linearization of the governing equations 
in the far field against the constant free-stream back- 
ground is relevant to many external flows.) Then, 
u = u — uo, where uq represents the corresponding 
free-stream parameters. Analogously to (1), we will 
assume that problem (2) is uniquely solvable and 
well-posed. 

Our final goal will be to find the solution u to (2a) 
on D^ ■ However, we cannot do it directly since 
system (2a) is not closed if considered separately. 
On the other hand, we also cannot directly solve 
on the computer neither (1) nor (2) because D ex is 
infinite. Consequently, we will need to supplement 
system (2a) by the special artificial boundary condi- 
tions (ABC’s) at T, so that the resulting problem 


= g(M), r ^An, 


<9u 

u, -X-, Vu,. 


rer, 


is closed and its solution on Di n is in a certain sense 
close to the solution of (2). To construct the ABC’s 
(3b), we will use the linearity of (2b)-(2c). 


1.2 Boundary Equations with Projections 

Clearly, the idea of exact ABC’s is that the re- 
lation (3b) should equivalently replace (2b)-(2c). 
Since the equations (2b)-(2c) are linear we can use 
the apparatus of generalized potentials and bound- 
ary projections in order to obtain the desirable 
boundary conditions. In so doing, Br from (3b) ap- 
pears to be some linear pseudo differential operator. 

The generalized potentials and boundary projec- 
tion operators that we employ for constructing the 
ABC’s were first proposed by Calderon 4 and then 
also studied by Seeley 5 . Ryaben’kii 7 (see also the 

O A 

recent paper and the book by Mikhlin, et al. ) 
had extended and modified the original approach, 
and also proposed an effective numerical technique 
for calculation of potentials and projections, this 
technique is called the difference potentials method 
(DPM). Additionally, Ryaben’kii in had for the 
first time shown how the generalized potentials and 
projections could be used for constructing the ABC’s 
for unsteady finite-difference problems. 

We, however, will for the moment restrict our- 
selves by considering the steady-state problems only, 
which means u = u(r) and u = u(r). To simplify the 
notations, we will further omit the tilde that distin- 
guishes between the “linear” and “nonlinear” solu- 
tions; it should not cause misunderstandings since 
hereafter we will mostly concentrate on the “linear 
part” of (2). 

On D in UD ex , we introduce the following auxiliary 
problem (AP) 


Lu = f , r € D in U D e 


suppf(r)c D in , (4a) 
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lu — .0, as r — > oo, (4b) 

which is the key element of our further construc- 
tion. In the theory of generalized potentials, the so- 
lution of the AP (4) plays approximately the same 
role as the convolution with the fundamental solu- 
tion plays in the classical potential theory. We will 
assume that problem (4) is uniquely solvable for any 
compactly supported right-hand side (RHS) f(r) and 
well-posed; the corresponding inverse (Green’s) op- 
erator is denoted G. For any function u(r), r € D ex , 
we also introduce its clear trace ^ on P, 

£r = Titu. (5) 


For the applications considered below (L in (2b) and 
(4a) is a second-order operator), £r is chosen as a 

u. — 

C is the normal to T. The corresponding space of 
vector-functions £r is called the space of clear traces. 
For any element of this space, we define the gen- 
eralized potential with the density 


Pfr= (g((Lw)| b J) 


D e 


(6) 


Here, w = w(r) is an arbitrary function with the 
only requirement of that it has the clear trace £p, 
Trpw = and the RHS for the AP, i.e., the func- 
tion that the Green operator G operates on in equa- 
tion (6), is given by 



on D in , 
on D ex . 


We also define the operator Pr, 


( 7 ) 


Pr£r = Tr r P£r, (8) 

as the trace (5) of the generalized potential (6); this 
operator obviously maps the space of clear traces 
onto itself. It turns out that Pp in (8) is a projection, 
Pp = Pp 5 it is called the Calderon boundary projec- 
tion. This operator will play a fundamental role in 
our further analysis. It is possible to show 6 ’ ' that 
those and only those £p that belong to the image of 
the projection, £p £ ImPp, i.e., satisfy the boundary 
equation with projection (BEP) 


Pr£r = £r, 


( 9 ) 


can be complemented from F to D ex so that the 
complement u — u(r), r £ D ex solves (2b)— (2c) . 
In other words, BEP (9) provides for a complete 
classification of those and only those densities £p of 
generalized potentials that are actually the traces 


of some solution u = u(r) to (2b)-(2c) on D ex . 
When the BEP (9) is satisfied, the corresponding 
complement on D ex can be restored as the potential 
(6). Clearly, since ImPp contains all those and only 
those £p that can be complemented on D ex so that 
the complement solves (2b)-(2c), then equation (9) 
equivalently replaces system (2b)-(2c), which means 
the BEP (9) serves as a desirable exact ABC of type 
(3b). We emphasize that these ABC’s are usually 
nonlocal; on the other hand, the algorithm for con- 
structing the ABC’s does not impose any essential 
limitations on the shape of T. 

BEP (9) provides for a general recipe to construct 
the exact ABC’s. However, any specific problem re- 
quires its own special approaches. First, certain as- 
sumptions in regard to the behavior of solution have 
to be done in order to pass from (1) to (2). For ex- 
ample, to linearize the governing equations one has 
to assume that the flow perturbations in the far field 
are small; this assumption can usually be verified a 
posteriori. Second, the AP (4), which needs to be 
actually solved for calculating Pp, is still formulated 
on the infinite domain. Special techniques must be 
employed to replace it by some finite-domain AP. 
Third, in practice the ABC’s are to be constructed 
for the discrete rather than for the continuous for- 
mulation of the problem, which means that some 
additional steps may be required for the actual im- 
plementation of BEP (9). Finally, the numerical 
process of the DPM ^ , which is used for the ac- 
tual calculation of ABC’s, may also vary in some de- 
tails from one problem to another. In Section 2, we 
will delineate the corresponding algorithm for the 
Navier-Stokes equations. Here, we will briefly de- 
scribe a simpler case of the two-dimensional external 
inviscid flows. 

1.3 An Inviscid Example 

Tsynkov 12 anc [ Sofronov and Tsynkov 13 con- 
sider a finite body (airfoil) immersed in an infinite 
flow of inviscid compressible fluid. The flow is gov- 
erned by the Euler equations and is assumed to be 
subsonic at infinity; for the purpose of numerical 
solution, the equations are discretized on a finite- 
difference O-type grid that is generated around the 
body. The computational domain Di n in 11 ’ 12, 13 
is formed by the grid; no special assumptions in re- 
gard to the shape of its external boundary T are 
done. Outside the computational domain, i.e., in 
D ex , the Euler equations are linearized against the 
free-stream background. Moreover, we assume the 
existence of the velocity potential in the far field 
and split the linearized system into elliptic (veloc- 
ity) and advection (entropy) parts. After the term 
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associated with the circulation of flow around the 
airfoil is subtracted, the regular part of the poten- 
tial of velocity perturbations can be described by the 
Prandtl-Glauert equation 


with zero boundary condition at infinity: 

<f(x, y ) — > 0, as x 2 + y 1 — ► oo. (10b) 


Equation (10a) can be easily reduced to the Laplace 
equation by means of an affine coordinate transform; 
in so doing, boundary condition (10b) obviously re- 
mains intact. To obtain the ABC’s for velocity com- 
ponents, we represent the solution to (10) in the 
form of a generalized potential and then construct 
the corresponding BEP. We use the AP formulated 
on a ring-shaped domain {Ri < r < i?o}; the exter- 
nal circle {r = Ro} encompasses P and the internal 
circle {r = Ri} lies inside Di n . At r — R\, we spec- 
ify homogeneous Dirichlet boundary conditions; at 
r = Rq, we specify boundary conditions 


d<f>k 

dr 


1 *| , 
+ — <Pk 
r=R 0 r 


r=R 0 


= 0 , 


k — 0, ±1, ±2, . . . , 


( 11 ) 


where <fk, k = 0, ±1, ±2, . . are the Fourier compo- 
nents of the potential (10) after the transform with 
respect to the polar angle. It is possible to make sure 
that this AP is uniquely solvable and well-posed for 
any R.HS concentrated inside the ring. Moreover, it 
is easy to see that once complemented to the exterior 
of the ring, the solution to this AP vanishes at infin- 
ity (compare to (4b)). Numerically, the AP is easy 
to solve by means of the discrete Fourier transform 
in polar coordinates. 

Analogously to the continuous formulations, the 
ABC’s in the discrete form should simply close the 
system that is solved inside Di n . For a second-order 
scheme employed inside the computational domain, 
we can always assume that the velocity components 
on the penultimate coordinate row P of the O-type 
grid are known, whereas the corresponding values on 
the outermost row Ti should be determined using 
the ABC’s. Therefore, we consider the penultimate 
coordinate line P as the actual artificial boundary. 
Referring the reader to the original work ’ ’ , 

as well as to Section 2, in which the Navier-Stokes 
case is analyzed, we skip here the details of the nu- 
merical procedure for calculating generalized poten- 
tials and projections and setting the inviscid ABC’s. 
We only mention that after the finite-difference BEP 


is obtained, it is solved so that the velocity on T pro- 
vided from inside Di n coincides with the gradient 
of the potential (10) in a certain generalized sense. 
Once the BEP is solved, we can find the discrete po- 
tential density for any data (velocity components) 
specified on T. When this grid density is known, we 
calculate the generalized potential and find the trace 
of its gradient on the outermost coordinate line Ti 
by means of interpolation; this procedure yields the 
ABC’s for velocities. Finally, the ABC’s for ther- 
modynamic parameters are obtained using local re- 
lations, specifically, the Bernoulli equation and the 
entropy advection equation. 

The technique of 1 10 for constructing the 

ABC’s was combined with an iterative method 14 by 
Sofronov for calculating steady solutions to the Euler 
equations. A few subsonic and transonic airfoil flows 
have been numerically studied; in the transonic case, 
local supercritical regions were always kept inside 
Di n so that one could treat the exterior linearized 
flow in D ex as purely subsonic. Numerical results 
presented in ^ demonstrate clear superiority of the 
nonlocal DPM-based ABC’s over the standard local 
techniques based on quasi-one-dimensional charac- 
teristic analysis. For a fixed computational domain, 
nonlocal ABC’s ' ’ provide for better accu- 

racy and faster convergence rate than local tech- 
niques; the entire numerical algorithm also appears 
to be more robust. Additionally, when the artifi- 
cial boundary approaches the airfoil, the solution 
obtained with the technique of ’ ’ appears 

to be essentially less influenced by the decrease of 
the size of computational domain than the solution 
obtained on the basis of the local boundary condi- 
tions. In other words, ABC’s 1 ’ allow one to 

maintain good accuracy of computations for much 
smaller computational domains than the standard 
boundary conditions do. These results, along with 
the geometric universality of ABC’s 1 ’ , make 

this approach useful for calculating external Euler 
flows. 

Analyzing this Euler example, we can see that the 
principle favorable circumstance that allowed us to 
relatively easy obtain the ABC’s is the availability of 
a finite-domain AP. This, in turn, is due to the fact 
that the Poisson equation admits variables separa- 
tion in polar coordinates. However, one obviously 
cannot rely on such circumstances in constructing 
the ABC’s for more general cases. For example, the 
linearized Navier-Stokes equations (see Section 2) 
and even the linearized Euler equations (without in- 
troducing the potential) most likely do not admit 
the separation of variables in any other coordinate 
system except in the Cartesian one. Therefore, we 
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need to develop some technique that would gener- 
ally allow us to approximate the solutions of the 
infinite-domain AP of type (4) by the solutions of 
a new finite-domain AP. In so doing, we should be 
able to control the corresponding approximation er- 
ror so that the resulting ABC’s be as close to the 
exact ones as desired. Below, the procedure for con- 
structing appropriate finite-domain approximations 
to the AP’s of type (4) is described for a particular 
class of systems of PDE’s with constant coefficients. 
For wider classes of systems, the possibility to use 
the same procedure is corroborated by the numerical 
experiments. 


where 8 > 0 does not depend on a. 

Let now U be a class of n-component vector- 
functions, in which we will be looking for the solu- 
tions u(a?, y) of system (12). The inclusion u(x, y) £ 
U holds for all those and only those vector-functions 
that satisfy the following conditions. 

• Each component of the vector- function u(#, y ) 
has continuous derivatives of all those orders 
that are involved in system (12). 

• Vector-function u(:c, y) can be represented as a 
Fourier integral: 


1.4 Systems of Simple Structure 

Hereafter in this subsection, we assume that all 
functions depend on only two space variables, i.e., 
r = (x, y), x and y being scalars. This assumption 
allows us to simplify the presentation and at the 
same time does not imply any loss of generality since 
the analysis of the case r = (x, y), y being a vector, 
is straightforward. 

Let us consider the system of PDE’s 

£ = A d?) u+f( * , ’' ) ' (12) 

— oo < x, y < oo, 


with respect to a n-component vector-function u = 
d 


u(x, y). In (12), A J is a n x n matrix, each 
entry of this matrix is a symbolic polynomial of 

d 

— — . The RHS f (a?, y) is a compactly supported n- 
oy 

component vector-function, f(®, y) = 0 for |ar — a\ > 
a, |y| > a, a being some positive constant. Sys- 
tem (12) can obviously be thought of as a particular 
class of systems of type (4a); consideration of only 
the first-order derivatives with respect to x presents 
no loss of generality since the higher derivatives can 
be eliminated by introducing additional variables. 
We designate by A (ia) the matrix that is ob- 


(0 


instead of 


tained by substituting ia into A ^ 

— . The entries of the matrix A(?’a) are, therefore, 

oy 

some polynomials of the real variable a (the coeffi- 
cients of these polynomials may, generally speaking, 
be complex). 


Definition 1 System (12) is called the system, of 
simple structure if the roots A? (a) of the equation 
det||A(ia) — AI|| = 0 (1 is the n x n identity ma- 
trix) satisfy the inequalities 


|3?Aj(of)| > 5(a) > 8 > 0, 


(13) 


DO 


u(z, y) = J u(x, a)e tay da , 

— OO 

(14a) 

oo 

«(*» a ) = J u (®> y) e iay dy. 

(14b) 


— OO 


• The derivatives involved in (12) can be repre- 
sented as 


dx 


CO 

h .1 [tC x ’ a \ 


0 iay 


da, 


(15a) 


OO 

^y) U = vfe J A ( ia ^ x ' a Y ayda • ( 15b ) 


• Fourier transformation u(ay a) of the function 
u(£, y) is bounded for each a, 

|u(ar, a)| < oo. (16) 

Note, |u(x, a) | in (16) can be thought of as the sum 
of the absolute values of components. 

Further, for any real Y > ‘2a we associate with 
(12) another system of PDE’s 

^r = A {^) UY+tYix ’ v) ' (n) 

— oo < x, y < oo, 


with the unknowns uy = uy(a?, y). In (17), fy (x, y) 
is a n-component vector-function that is periodic in 
the y direction with the period Y and that coincides 
with the RHS f(x, y) of system (12) for |y| < Yf 2. 

Let Uy be a class of n-component vector- 
functions, in which we will be looking for the so- 
lutions uy(£, y) of system (17). The function 
uy(x, y) belongs to Uy if and only if it satisfies the 
following conditions. 
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• Each component of uy(®, y ) has continuous 
derivatives of all those orders that are involved 
in system (17); moreover, uy((E, y ± Y) = 
uy(s, y). 

♦ Vector-function uy(a?, y) can be represented as 
a Fourier series: 


k—oo 


uy(x, y) = 


uy^) 


1_ 

Y 


k — — oo 

(18a) 

Y/2 

/» 


■ / uy(z, y)e~ lky ^dy. 

(18b) 

-Y/2 



• The derivatives involved in (17) can be repre- 
sented as 


chi y 

dx 


k—oo r 


= E 


k — — oo 


TC Y - (x) 


Jkyl? 


A1 %, |Uy = 


(19a) 

(19b) 


k = oo 


E A 


k = — oo 


. 27 Tk \ 


7 — j n Yk {x)e lky y . 


• Fourier coefficients uy fc (a?) of the function 
uy(#, y) are bounded for all k, 


|uy fc («)| < oo. 


(20) 


Theorem 1 The system of simple structure (12) 
may have at most one solution u(a?, y) E U , and sys- 
tem (17) may have at most one solution uy(r, y) E 

U Y . 


Proof. Clearly, it is sufficient to show that the 
homogeneous system 


<9u 

dx 


= A 



u 


has only trivial solution u(x, y) = 0 in the class 
U . Let u(x, y) E U be some solution to this ho- 
mogeneous system. We represent this solution as a 
Fourier integral in accordance with (14a) and, using 
(15), obtain 



— u(x, a) — A(ior)u(a:, cr) 
dx 


e iay da 


0. 


Therefore, for any a 


du(x, a) 
dx 


= A(ia)u(x, a), 


( 21 ) 


and u(a?, a) is bounded. On the other hand, all roots 
of the characteristic equation of system (21) always 
have non-zero real parts since (12) is a system of 
simple structure. It is well-known that in this case 
system (21) has a unique bounded solution u(a:, a), 
which is identically zero. Consequently, 


u(a?, y) = 



/ u(x ’ a)e ‘ ayda s 0 


— OO 


The statement about uy(x, y) can be proven anal- 
ogously. □ 

To formulate the next theorem, we intro- 
duce the following notations. Let <p(x, y) be 
some n-component vector- function, <p(x, y) = 
(<P!(x, y), <p 2 (x, y), • • • , <pn{x, y)) T . Designate 


|y>(*. y) I = M*’ y)\> 

j - 1 


d p+q <p(x, y) 
dx p dy q 


n 


E 


d p+g (fj(x, y) 
dx p dy q 


IK*, y)IU 


E 

P+<7=0 


sup 

y 


d p+q (p(x, y) 
dxPdy q 


(22a) 

(22b) 


(22c) 


Theorem 2 Consider an arbitrary bounded domain 
D C R 2 . Let z be an arbitrary non-negative integer 
number. Then, one always can find a sufficiently 
large number k — k(z, D ) such that for any com- 
pactly supported function f(x, y), f(x, y) = 0 for 
\x — a\ > a and |y| > a, that satisfies ||f [j^ < oo 
(see (22)) the solutions uEf and uy E Uy of sys- 
tems (12) and (17), respectively, exist (Y > 2 a is 
arbitrary ) and 


II u - u HL,d < 


const 

Y K ’ 


(23) 


where k — k(z, k , D), and k 


oo as 


The importance of Theorem 2 is that it allows us 
to approximate the solution of the infinite- domain 
AP by the solutions of another problem, which is 
periodic and, therefore, finite, in all but one Carte- 
sian directions. Recall, to calculate the ABC’s we 
need to know the projection Pp (see (8)), which is 
the superposition of the trace Trp (see (6)) and the 
potential P (see (5)). We, therefore, conclude that 
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in order to obtain the ABC’s we do not need to know 
the potential, i.e., the solution of the AP, everywhere 
on D ex \ we need to know it only on some neighbor- 
hood of T. Consequently, the approximation of u 
by u y on any finite domain D as the period Y in- 
creases (see (23)) is sufficient for achieving our goal 
of constructing the ABC’s. 

The proof of Theorem 2 is essentially based on 
the estimates for one-dimensional fundamental so- 
lutions obtained by Godunov and Gordienko in 
Specifically, let 

^ = Av + f(D) (24) 

be an abstract system of n ordinary differential equa- 
tions (ODE’s) with constant coefficients, and f (x) be 
a continuous bounded function. Let also the roots 
A of the characteristic equation det|A — AI| =: 0 of 
system (24) satisfy the inequality 

|3t\| > 8 > 0. (25) 

We supplement system (24) by the condi- 
tion of boundedness of its solution v(a)) = 

(^i(x), V 2 (x), . . . , v n (x)) on the entire line — oo < 
x < oo, 

1 1 v(m) | J < const. (26) 

Because of the equivalence of norms on a finite- 
dimensional space for a fixed n, the choice of the 
specific norm ||v(a?)|| in (26) is not essential. Here- 
after, we will always use the l\ norm of vectors, i.e., 
consider as a norm of a vector the sum of absolute 
values of its components. We now formulate the fol- 
lowing result, see Fundamental solution G(x) of 
the corresponding linear differential operator from 
(24) exists and is unique in the class of functions 
that meet boundary condition (26). This fundamen- 
tal solution is actually a n x n matrix that satisfy 
the inequality 

( II A II \ n ~l 

U) e~ 6 ^, (27) 

where f2(n) is some number that depends only on 
n, and the norms of the matrices G(x) and A are 
consistent with the chosen vector norm. Note, ac- 
cording to the definition of the fundamental solution, 
one can represent the solution to (24), (26) as 

OO 

v(a?) = f G(x — t)f(t)dt. (28) 


We will be looking for the solution of system (12) 
in the form of the Fourier integral (14a), where the 
Fourier transformation u(a?,a) is subject to the de- 
termination. Analogously to (14), we represent the 
RHS f(x, y ) of system (12) as 

OO 

i(x, y) — —j= f f(x, a)e iay da, (29a) 

v2tt J 

— OO 
OO 

f(«> <*) = J f(s, y)e~ iay dy. (29b) 

— OO 

Substituting (14a) and (29a) into (12) and using 

(15) , we formally obtain for u(x,a) the following 
system of ODE’s with constant coefficients 

du(ar, n) .. . „ . . ~ , „ 

— — — - = A(?o;)u(^, a) -j- f(x, a). (30) 

System (30) depends on a as on the parameter. We 
will supplement this system by boundary condition 

(16) . 

To make sure that the function u(x, y) obtained 
using (14a) (provided that u(«,a:) satisfies (30), 
(16)) does belong to U and does solve (12), we will 
first have to make certain assumptions in regard to 
the smoothness of the RHS f(x, y) (recall, this func- 
tion is always supposed to have compact support). 
Then, we estimate from above the Fourier transfor- 
mation f(x, a) (see (29a)) and its derivatives by cer- 
tain rational functions of a. Using the estimates for 
f(x,a), inequality (27), and representation (28), we 
can estimate from above the solution u(x, a) of (30), 
(16) and its derivatives with respect to x by the (22)- 
type norm of f(;r, y) multiplied by a certain rational 
function of a and a certain exponential function of x. 
Further, the differentiability of u(^,a) with respect 
to a and the estimates for the derivatives (in partic- 
ular, the rate of decay for big |a|) can be obtained 
by induction with respect to the order of differenti- 
ation. Finally, using the aforementioned estimates 
we make sure that the inverse Fourier transforma- 
tion (14a) does exist in U and therefore, provides 
for the solution to (12); using the same estimates, 
we also make sure that inequality (23) does hold, 
first for z = 0 and then for z > 0. 

This brief outline will be transformed into the full 
proof of Theorem 2 in a future paper. Here, we will 
show how to handle the last remaining infinite direc- 
tion, x , and to therefore finally replace the infinite- 
domain AP by a new problem formulated on some 
bounded domain. 
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Since we have already replaced the original 
infinite-domain formulation of the AP by the peri- 
odic formulation with respect to y , then, instead of 
solving systems (30), (16) for each a, — oo < a < oo, 
we will need to solve the systems 

Y k ( x ) = A uy fc (z) + fy,(:c) (31) 


with boundary condition (20) for each k, k = 
0,d=l,±2, ... Note, system (31) is obtained on the 
basis of representation 


f y(x,y)= 

k = — oo 
Y/2 

*Y k (x)=y f f Y (x, y)e~ tky ^dy 
-Y/2 


(32a) 


(32b) 


by substituting (18a) and (32a) into (17) and using 

( 19 ). 

For each k, k = 0,±1,±2,..., system (31) for- 
mally needs to be solved on the entire line — oo < 
x < oo. However, since we need to know the solution 
of the AP only on some finite neighborhood of T, we 
can always choose a sufficiently large but still finite 
segment on the x axis, on which it would be suffi- 
cient to calculate the solution to (31). Without loss 
of generality, we may think that it is the segment 
(0, X), A' > 2a. Since the RHS f(a?, y) is compactly 
supported, system (31) is homogeneous for x < 0 
and x > X . To ensure that boundary condition (20) 
is met, we have to prohibit for x > X all the eigenso- 
lutions of the homogeneous counterpart to (31) that 
increase as x — > +oo and to prohibit for x < 0 all 
the eigensolutions that increase as x — » — oo. This 
can be done by imposing the boundary conditions 


n 

3£AjO)>0 


(A* - a .,•(*)!) 


uy„(0) = o 


and 


(33a) 


n (a* 

dtXj(k)<0 


A j(fc)I) 


un(A') = o 


(33b) 


uy fc (a?) vanishes at infinity (which is a stronger prop- 
erty than simply boundedness (20)). 

Thus, we have shown how to replace the infinite- 
domain AP for a system of simple structure by the 
new problem formulated on the rectangle (0, X) x 
(— y/2, y/2) for the same RHS so that the solutions 
of these two problems are arbitrarily close to one an- 
other on a fixed finite neighborhood of Di n . We only 
have to mention, that the “asymmetric treatment” 
of the Cartesian directions x and y is caused mostly 
by the physical concerns. Specifically, typical ex- 
ternal flow formulations (see Section 2) would bear 
some natural non-isotropy when we treat the Carte- 
sian direction x as a streamwise and the Cartesian 
direction y as a cross-stream. In so doing, exact 
analytic treatment (33) of the Cartesian direction 
x seems to be most relevant for practical computa- 
tions, since the periodic treatment for x may result 
in much larger values of period required for achiev- 
ing the same accuracy than the periodic treatment 
for y. On the other hand, for the purposes other 
than the practical computation of projections, we 
may consider Fourier representation of the solution 
to the AP in all Cartesian directions. In particular, 
this has been done by Tsynkov in when proving 
the solvability of the AP for the linearized thin-layer 
equations in the sense of tempered distributions. 

To conclude the introductory part of the paper, 
we will also briefly comment on the relation between 
the systems of simple structure and other equations 
that can often be encountered in practice. 


1.5 Wider Classes of Systems 

Many systems of equations that are frequently 
solved in mathematical physics may not appear to 
be the systems of simple structure. For example, 
consider the system 


du 

dv 

dx 

y G R , 


d 2 ii 

— +/m + /0 e, y), 

y — const, 


(34) 


which is obtained when replacing the equation 

£ + 0-/- = /(*•») W 


at x = 0 and x — A, respectively. In (33), Xj(k) 
are the eigenvalues of A k = A and the ma- 

trix products are computed taking into account the 
multiplicities of these eigenvalues. Clearly, since for 
the system of simple structure |3£Aj > 6 > 0, bound- 
ary conditions (33) actually imply that the solution 


by the first-order system with respect to x. For (34), 
the matrix A(ia) becomes 


0 1 
a 2 + y 0 ’ 

and its eigenvalues are A(») = ±\/ a' 2 + y . When 
y > 0 we have |3ftA(a)[ > ^/y > 0, and system 


A(io;) = 
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(34) is a system of simple structure. However, when 
y = 0 (in this case, equation (35) is the Poisson equa- 
tion) or when y < 0 (equation (35) is the Helmholtz 
equation), the real part of A(a) for some cv’s is zero, 
which means that system (34) is not a system of sim- 
ple structure. Accordingly, the Poisson equation on 
R? not always has a bounded solution. As for the 
Helmholtz equation, the condition of vanishing of its 
solution as r = (# 2 + y 2 ) 1 / 2 — *■ oo does not ensure 
the uniqueness. To select the unique solution, one 
has to impose a more fine Sommerfeld radiation con- 
dition, which, in particular, also implies vanishing of 
the solution at infinity. 

As shown above, if the linear differential opera- 
tor L from (2b) being represented in the form of 
(12) implies the inequality (13), i.e., the correspond- 
ing system of PDE’s meets the Definition 1 of the 
systems of simple structure, and if the boundary 
conditions at infinity (2c) are actually the condi- 
tions of a sufficiently fast decay, then the ABC’s 
at the artificial boundary T = dDi n can be con- 
structed using the finite-domain AP on the rectan- 
gle (0, X ) x (— "T/2, Y j 2). Indeed, under these con- 
ditions the infinite-domain AP (4) formulated for 
the RHS (7) can be replaced by the finite-domain 
AP using Theorem 2 and boundary conditions (33). 
If, however, system (2b) represented in the form of 
(12) is not a system of simple structure, then we 
can sometimes approximate its solutions by the so- 
lutions of a specially chosen system of simple struc- 
ture. (Note, the representation in the form of (12), 
if possible at all, may require the solution with re- 
spect to dw/dx after the Fourier transform; in so do- 
ing, the entries of A(«o;) may become rational rather 
than polynomial functions of a). 

Let B = B — , £ ) be a n x n matrix, which ap- 

\dy J 

proaches the zero matrix as e — *• +0, and let the 
roots A(a, e) of the equation 

det||A(*a) + B (ia, e) — A(a, 6r)I| j = 0 (36) 

satisfy for each e the inequality 
\m(a, e)\ > 8(e) > 0. (37) 

Let, in addition, the problem 


$(u, Vu,...,£) = g(i,i/), (x, y) e D in (38a) 


du(x, y, e) 
dx 

(x, y) G D ei 


u(ar, y, e) — *• 0, as r — *• oo 


(38c) 


be uniquely solvable and let its solution u(a?, y, e) 
approach the solution to the steady-state counter- 
part of problem (2) as e — *• +0 on any finite sub- 
domain of R 2 . Then, for constructing the ABC’s on 
T one obviously can use the AP 


du(x, y, e) _ 
dx 

A (i) +B (4 >£ )] u+,( * ,I ' ) ’ 

(x, y) € R 2 , suppf(®, y) C D in 


(39a) 


u(x, y, e) — *• 0, as r — *• oo, (39b) 

where e is chosen sufficiently small. 

For example, if y — 0 in system (34), then the 
correction 




will transform (34) into a system of simple struc- 
ture. In the case of the Helmholtz equation (y < 0 
in (34)) supplemented by the Sommerfeld radiation 
conditions at infinity, we choose 



0 0 ' 
is 0 


This choice corresponds to the well-known principle 
of limitary absorption, which allows one to approx- 
imate both the Helmholtz equation and the Som- 
merfeld radiation conditions with the increasing ac- 
curacy as e — ► +0. We should note that another 
choice of B, 



0 0 ‘ 
—ie 0 


also provides for the system of simple structure 
(39a). However, in this case we approximate (as 
e — — * +0) the solution composed of incoming waves 
rather than the one that meets the radiation condi- 
tions. 

Finally, we mention that for a given system of 
PDE’s it is, generally speaking, not always easy to 
find a good analogue to the principle of limitary ab- 
sorption so that the solutions of this system can be 
approximated by the solutions of a certain system 
of simple structure. However, the replacement of an 
infinite-domain AP by the periodic problem may still 
be possible, the validity of such a replacement can 
be verified a posteriori by means of the numerical 
experiments as done, e.g., in our work . 
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2 ABC’s for External Viscous Flows and 


2.1 Governing Equations 

We consider a flow of compressible viscous fluid 
over a finite body or configuration of bodies (e.g., 
single-element or multi-element airfoil in two space 
dimensions). The flow is supposed to be uniform and 
subsonic at infinity. As mentioned above, we intend 
to actually calculate the flow only on some finite 
computational domain D{ n (which contains the im- 
mersed body(ies)). The infinite exterior D ex of the 
domain Di n is truncated, and the equations there 
along with the boundary conditions at infinity are 
to be replaced by the ABC’s at T = dDi n . 

To construct the ABC’s, we linearize the gov- 
erning equations in D ex against the constant free- 
stream background. Generally, to justify the pos- 
sibility of linearization we have to make sure that 
first, the flow perturbations in the far field are small, 
and second, the differential equations that describe 
these small perturbations' are linear. As concerns the 
first assumption, it does hold for the external viscous 
flows, at least when the size of Di n is much larger 
than the size of the immersed body. The linearity 
of the governing equations for small perturbations 
is easy to establish for the subsonic flows, i.e. , when 
the free-stream Mach number Mo is not too close to 
one. For the transonic flows, it basically requires a 
special study (see our work 18 and also the recent 
paper 1 ^). The situation appears different for three 
space dimensions, when the far field is essentially lin- 
ear, and for two space dimensions, when the consid- 
eration of a simplified potential model may formally 
require to introduce some nonlinear corrections in 
the transonic limit, i.e., as Mo — > 1. We, however, 
always consider the full flow system rather than the 
potential equation, we also never approach the tran- 
sonic limit too closely. In so doing, we use the far- 
filed linearization for two space dimensions as well. 
In practice, the validity of the far-field linearization 
is always verified by an a posteriori numerical check. 

Either one of the two following systems of PDE’s 


dp du dv _ 

dx "I" dx dy 

du dp 1 '4 d 2 u 

dx dx Re 3 dx 2 


dv dp 1 '4d 2 v 

dx ^ dy Re _3 dy 2 
dp 1 dp 7 
dx Mq dx Re Pr 


1 d 2 v d 2 u 

3 dxdy dy 2 _ 

1 d 2 u d 2 v' 

3 dxdy + dx 2 _ 
d 2 p d 2 p 

JhP + dtf~ 


(40) 


= 0 
= 0 


1 / d 2 p d 2 p\ 

7 M 2 \ dx 2 + dy 2 ) 


= 0 


dp du dv 
dx^dx^ dy 



du dp 1 d 2 u 

dx dx Re dy 2 


dv dp 14 d 2 v 

dx ^ dy Re 3 dy 2 

0 


dp 1 dp 7 

'd 2 p 

1 d 2 p 

dx Mq dx Re Pr 

.dy 2 

7 Mq dy 2 


can be used (and has actually been used) for con- 
structing the ABC’s in two space dimensions for 
steady-state flows. Hereafter, we will concentrate 
mostly on this specific case (2D steady-state flows); 
in the end of the paper, we will also present some re- 
cent three-dimensional results (see our work 
as well as the earlier papers ^) and comment 
on the possible approaches to treating the time- 
dependent problems. 

System (40) represents the linearized dimension- 
less full Navier-Stokes equations, and system (41) 
represents the linearized dimensionless thin-layer 
equations. In (40) and (41), u, v, p , and p are the 
perturbations of the Cartesian velocity components, 
pressure, and density with respect to the correspond- 
ing free-stream parameters, Mo is the free-stream 
Mach number, Re is the Reynolds number, Pr is 
the Prandtl number, and 7 is the ratio of specific 
heats. In both cases we assume that the direction 
of flow at infinity coincides with the positive x di- 
rection, and that the gas is perfect. To obtain the 
dimensionless quantities, we use the following scales: 
uq, for velocity components; p 0 , for density; poUo 2 , 
for pressure; po, for viscosity; characteristic size L 
(typically, airfoil chord), for all distances. Here, the 
subscript “ 0 ” denotes the free-stream parameters. 

As mentioned above, to construct the ABC’s we 
need to be able to solve the AP for a nonhomoge- 
neous counterpart to either (40) or (41) driven by a 
certain compactly supported RHS. In our work 16 
we have, in particular, shown that when supple- 
mented by a compactly supported RHS, system (41) 
is always solvable in the sense of tempered distri- 
butions (see, e.g., book 22 by Hormander or 28 by 
Vladimirov for a detailed description of the concept 
of distributions and its application to solving the 
PDE’s). We have also shown that if this solution 
satisfies the boundary condition 

{u, v , p, p) — r ( 0 , 0, 0, 0 ), 

as r = {x 2 + y 2 ) 1 / 2 ► + 00 , 


( 42 ) 
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then it is unique in the class of distributions vanish- 
ing at infinity. Note, the solvability in S' (space of 
tempered distributions) does not necessarily mean 
the fulfillment of (42). To make sure condition (42) 
does hold, we need to do certain assumptions about 
the RHS £(x, y ). From the very beginning, we al- 
ways think that f is compactly supported and that 
f G L 1 (R 2 ), which is no loss of generality. If we 
additionally assume that f G L 2 (R 2 ) (which basi- 
cally presents no loss of generality as well), then the 
solution u = (u, v , p, p) to the infinite-domain AP 
(of type (4)) can be represented asu = u^'+ u( 2 ), 
where u 1 ^ 1 ) satisfies (42) and G L 2 (R 2 ); the lat- 
ter inclusion can be treated as “a generalized decay 
at infinity”. If, however, we require that f(;r, y ) be 
sufficiently smooth on R? so that its Fourier trans- 
formation with respect to both x and y belongs to 
L 1 (R 2 ), then we can show (see ^ 6 ) that the solu- 
tion u to the AP for system (41) meets boundary 
condition (42). Similar results for three space di- 
mensions have been obtained in our work iy . Let 
us also note that the work 16 is devoted to studying 
a wider class of formulations than only the steady- 
state flows, specifically, we study there the flows that 
oscillate in time. The aforementioned solvability re- 
sults for the steady-state thin-layer equations can 
be obtained as one consequence of the considerations 
of lf \ We will briefly review the general results of ^ 
in the end of this paper. 

2.2 Geometric Setup 

Let us now introduce the geometric setup typical 
for external flow problems in two space dimensions; 
an example is shown in Figure 1. We are interested 
in calculating the flow around an airfoil, the single- 
element configuration presented on the figure does 
not imply any loss of generality since the treatment 
of external boundary for the multiple immersed bod- 
ies would basically be the same. To calculate the 
flow, we first generate the grid around the airfoil; 
for this specific example it will be a C-type curvi- 
linear grid, which actually forms the computational 
domain Di n . The nonlinear flow equations are dis- 
cretized and solved on this grid inside D 2n . How- 
ever, analogously to the continuous case (see Sec- 
tion 1) the discrete system inside D ?n is subdefinite 
unless we supplement it by some ABC’s. Indeed, 
the stencil of the finite-difference operator used in- 
side Di n cannot, generally speaking, be applied to 
those nodes of the C-grid that are located near the 
external boundary, e.g., it cannot be applied to any 
node of the outermost coordinate row of this grid, 
since in so doing the part of the stencil may sim- 
ply “fall out” of the domain. Therefore, the discrete 


system inside D in without ABC’s would merely have 
less equations than it has unknowns. Consequently, 
unlike the continuous case, for which to close the 
system inside Di n means to set the ABC’s exactly 
at the continuous external boundary, to close the 
system inside Di n in the discrete framework means 
to provide for some additional relations between the 
values of the solution in the nodes located in a cer- 
tain external part of the grid. For example, if the 
scheme employed inside Di n is written on the 3x3 
stencil, which, in particular, corresponds to a widely 
used second-order central-difference approach, then 
the ABC’s should provide for the missing relations 
between the values of the solution on the penulti- 
mate and outermost coordinate rows of the C-grid. 
On Figure 1, the penultimate coordinate row is des- 
ignated T and the outermost row is designated Fi. 

Henceforth, we will treat the penultimate coor- 
dinate row T of the C-grid as a formal continuous 
artificial boundary. This, in particular, means that 
the outermost curve Ti belongs already to the area, 
in which we linearize the governing equations. It also 
means that if we were looking for the continuous so- 
lution on Di n , then we would need to construct the 
ABC’s exactly at the penultimate coordinate line L. 
For the discrete formulation, however, we need to 
obtain missing relations between the values of the 
solution on T and Ti (see above). To do that, we 
will first formulate the ABC’s of type (9) exactly at 
T and then use the generalized potential (6) for com- 
plementing the boundary data from T to D ex , the 
trace of this complement on T i will provide us with 
the unknown values of the solution at the outermost 
coordinate line of the C-grid. In so doing, we not 
only obtain the desirable missing relations that close 
the discrete system in Di n , but at the same time au- 
tomatically make sure that these relations are right 
in the sense that they properly take into account the 
structure of the solution in the far field; the latter is 
true because the complement we construct on D ex 
is obtained on the basis of (9). 

2.3 Computation of the Potentials and 
Projections 

In fact, we, of course, cannot calculate directly the 
continuous generalized potentials (6) and boundary 
projections (8); instead, we calculate their discrete 
counterparts called the difference potentials and the 
difference boundary projections, respectively. The 
corresponding numerical procedure is based on ap- 
plication of the DPM ^ . The issues of consistency 
and convergence for the difference potentials and 
their continuous prototypes have been studied by 
Ryaben’kii in and Reznik in . 
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Figure 1. Configuration of domains. 


As described in Section 1, the AP for the nonho- 
mogeneous version of either system (40) or (41) is 
first formulated on the entire plane and then trun- 
cated so that one needs to solve it only on the rectan- 
gular domain Dy — (0, X) x (— Y/2, Yj 2); this aux- 
iliary domain should fully contain both T and Fi , see 
Figure 1. We will now describe the finite-difference 
formulation of the AP and the DPM-based algo- 
rithm for calculating the difference potentials and 
projections. 

Let us introduce in Dy two Cartesian grids, A4° 
and Af°. The grid M° will be used for specifying 
the RHS for the finite-difference AP, the grid J\f° is 


the one, on which the solution to this AP will be 
defined. We also introduce the space F° 3 f° = fjvjo 
of the RHS’s for the AP, the space U° 3 u° = Ujyo 
of its solutions, and the finite-difference operator 
hh ■ U° i — > F°, which can be a discretization of 
the left-hand side of either (40) or (41). Note, in 
the work ’ Z0, ZD we have used the operator L* 
obtained by the second-order central-difference ap- 
proximation of (40), in so doing the grids A4° and 
Af° coincide except on the lines x = 0 and x = X , 
where the RHS grid A4° is simply not defined. In a 
later work ’ ’ , we have used the operator L h 

obtained by the second-order approximation of (41) 
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with the central differences along y and the first- 
order differences along x. In so doing, the RHS grid 
M° is shifted in the x direction with respect to Af° 
by the half of the grid size, and the finite-difference 
equations of the AP are written as 

L t u° = D ^ + hLZ^ hL + (43) 

tlx 


Ip 1 ( Um >-7 + 1 u m+l,j + l u m+l,j-l i , 

2 V 2 hy 2hy ' 


1 jj ( u mj+l 2u m,j T U m j' _ i 


h l 

2llm+l.? T U 


T 




h? 

U y 


= f, 


m+l/2,j > 


m = 0, . . . , M — 1, j — 0, . . . , 2J, 


where h x and h y are the Cartesian grid sizes; sub- 
script m corresponds to the Cartesian direction x, 
and M - f 1 is the total number of nodes of the grid 
Af° in this direction; subscript j corresponds to the 
Cartesian direction y, and 2J + 2 is the total number 
of nodes of the grid JC° (and A4°) in this direction. 
The 4x4 matrices D, F, and H are given by 


D = 


10 0 1 

10 1 0 

0 10 0 

0 0 1 -M 0 “ 2 


0 10 0 
0 0 0 0 
0 0 10’ 
0 0 0 0 


(44) 


H = - 


Re 


0 0 0 0 

10 0 0 

04/30 0 

0 0 yPr" 1 -Pr~ l M 0 ~ 2 


Note, we omit the superscript “0” for the solutions 
and the RHS’s of the AP when referring these func- 
tions to specific grid nodes, see (43). 

Following the considerations of Section 1, we re- 
quire that any function u° £ U° be periodic in the 
y direction with the period Y, i.e., 

u m , 0 = u m>2 j+i, m = 0,...,M, (45) 

u m _i=u mi2 j, m = 0, 


As concerns the RHS’s f° £ F°, they may, generally 
speaking, differ from zero only for those nodes of 
the grid A4° that belong to Di n . Instead of (18) 


and (32), we now use the standard discrete Fourier 
transform (see, e.g., 16, 17, 25 for details), and obtain 
the following family of systems of ordinary difference 
equations 

= fm+l/2,fc ■> (46) 


m = 0, . . . , M — 1, k = — J, . . . , J, 

which represent the finite- difference analogue of 
equations (31). In system (46), subscript k is the 
dual discrete Fourier variable that corresponds to 
the Cartesian direction y, and f m +i/2,/ t are 

the discrete Fourier transformations of the solution 
and the RHS, respectively; A* and are the square 
coefficient matrices, which depend on k but do not 
depend on m. These matrices are obviously deter- 
mined by the structure of the finite-difference op- 
erator L/j. For system (43), Ak and B& are given 
by 

A * = i D + f F + t H ’ (47 > 


„ 1 Tif tk 

B ‘ = -AT D+ y r+ ¥ H ' 


where the matrices D, F, and H are defined in 
formula (44), = isin (kh y ^- \ /h y , and tk = 



Specific expressions for A& 


and Bfc that correspond to system (40) can be found 
in and those that correspond to a more general 
case of time-periodic flows can be found in 16 . Note, 
the matrices A* and B^ from (47) have order 4, 
whereas for system (40) these matrices would have 
order 8 (see 17 ), this is caused by the necessity to 
introduce additional variables when reducing the or- 
der of differencing with respect to x from the second 
to the first. 

To complete the formulation of the finite- 
difference AP on the rectangle Dy , we have to spec- 
ify the boundary conditions at x = 0 and x = X, 
i.e., at m — 0 and m — M . Analogously to how it 
is done in Section 1 for the continuous formulation, 
we first note that system (46) formally considered 
on the infinite mesh — oo < m < oo is homogeneous 
for m < 0 and m > M . Then, to meet boundary 
condition (42), we have to prohibit all the growing 
eigensolutions of (46) for m — » — oo, as well as for 
m — *■ Too. This can be done by imposing for each 
k, k = — J, . . . , J , the boundary conditions 


n 

|Ar(*)|>l 


(Qi-A r (fc)I) 


Uo ,fc = o, 


(48a) 
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n 

|A r O)!<i 


(Q* - K(k)l) 


UM,fe — 0, 


(48b) 


where Qk — (A *,) -1 B*,, A,.(&) are the eigenvalues of 
Qfc, I is the identity matrix, and the matrix prod- 
ucts are calculated in accordance with the multi- 
plicities of eigenvalues. Discrete boundary condi- 
tions (48) are analogous to the continuous boundary 
conditions (33). The only difference is that unlike 
(33a.) and (33b), equations (48a) and (48b) are not 
symmetric. In formula (48b) (downstream bound- 
ary condition), we admit the eigenvalues |A(&)| = 1, 
which in the continuous case would correspond to 
|9?A(a)| = 0. The reason for this asymmetry is 
that none of the systems (40) and (41) is actually 
a system of simple structure. For a — 0 (k — 0 in 
the discrete formulation), we have (multiple) eigen- 
value with |3?A((*)| = 0 (|A(fc)| = 1 ) for both (40) 
and (41). Unfortunately, neither for the linearized 
Navier-Stokes nor for the linearized thin-layer equa- 
tions we are unaware of any good analogues of the 
principle of limitary absorption (see Section 1) that 
could have helped us to approximate the solutions 
of these systems by the solutions of systems of sim- 
ple structure. We have, therefore, to entirely rely on 
the numerical experiments as the means to justify 
the possibility to replace the original AP by the pe- 
riodic one. The computations that do corroborate 
the possibility to introduce the periodic formulation 
are reported in . Moreover, it is also possible to 
make sure that the matrix Qk for |A(/?)| = 1 (k — 0) 
still has a full system of eigenvectors, which means 
that the corresponding eigensolutions of the homo- 
geneous one-dimensional system are at most oscilla- 
tory, but never increasing. Using this fact, one can 
show (see ^) that the resulting solution obtained by 
means of an inverse Fourier transform still vanishes 
at infinity, even if for a selected finite number of a’s 
(fc’s) we require only the boundedness (48b) rather 
than the true decay of the one-dimensional solution 
in the Fourier space. 

Thus, we have finally completed the formulation of 
the finite-difference AP, and have, therefore, defined 
its Green (i.e. , inverse) operator G/ t , Gh : F° \ — > 
U°. It is easy to see that this AP is uniquely solvable 
for any compactly supported RHS and well-posed, 
the well-posedness can be established on the basis 
of the considerations of 29 . An effective numeri- 
cal algorithm for solving the AP (more precisely, for 
solving one-dimensional systems (46) with boundary 

on 

conditions (48)) is described in our work 130 . This 
algorithm can be referred to as a version of the well- 
known successive substitution technique, but with- 


out its “inverse” or “resolving stage”. The partic- 

on 

ular efficacy of the approach of is based on the 
fact that boundary conditions (48) are formulated in 
terms of eigen subspaces of the operator Q&. 

Let us now split the nodes of the grid A 4° into two 
groups, Min = and M ex = 

Then, we apply the stencil of the finite- difference 
operator L/j to each node of the set Mi n and call 
the union of all these stencils A*„. Analogously, we 
apply the stencil of L h to each node of M ex and ob- 
tain Af ex ■ The intersection of the two sets and 
A f ex is called the grid boundary 7 , 7 = Am C\Af ex - 
The subset 7 C Af° is actually a multi-layered fringe 
composed of those nodes of the grid Af° that are lo- 
cated in a certain sense near the continuous artificial 
boundary F. 

For any function u° £ U° we introduce its differ- 
ence clear trace £ 7 on the grid boundary 7 as merely 
a contraction, 

£7 = Tr 7 u 0 = u^o| r (49) 


Since 7 is a multi-layered set of nodes, £ 7 of (49) in 
a certain sense models of (5). Let us now intro- 
duce the space of difference clear traces that would 
contain all grid vector-functions of dimension 4 de- 
fined on 7 . For each element £ 7 of this space, we can 
construct the generalized difference potential 


P«tr= (^(WL,.)) 



(50) 


where w° £ U° is chosen so that it has the 
trace £ 7 , Tr 7 w° = £ 7 , and is arbitrary in the 
rest. For example, one always may choose w° = 

{ c on 

7 ’ . rn , . As concerns the RHS for the 

0 , on A u \7 

difference AP, i.e., the function that the discrete 
Green operator G^ operates on in equation (50), it 
is defined as 


( L * w °)L, 


LftW 0 , on Min , 
0 , on M ex . 


(51) 


Clearly, to calculate the difference potential P ea; £ 7 
(50), which is analogous to the continuous gener- 
alized potential ( 6 ), we need to actually solve the 
difference AP. 

Finally, we define the operator P 7 , 


P t £ 7 = Tr 7 P ea; G 


(52) 


as the composition of potential (50) and trace (49). 
This operator obviously maps the space of difference 
clear traces £ 7 onto itself. As in the continuous case 
(compare (52) with ( 8 )), P 7 appears to be a projec- 
tion, P^ = P 7 , it is called the difference boundary 
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. *7 

projection. It is possible to show D| ' that those and 
only those £ 7 that satisfy the BEP 

P 7 £t=«t, (53) 

i.e., belong to the image of the boundary projec- 
tion P 7 , £ 7 E ImP 7 , can be complemented on Af ex 
so that the complement UjV ex satisfies the homo- 
geneous equation L^uy ei — 0 and boundary 
conditions (45), (48) of the difference AP. In other 
words, BEP (53) provides for an exhaustive classifi- 
cation of those and only those £ 7 ’s that can be repre- 
sented as a trace of some solution ujy ex of the equa- 
tion L/j u = ^Mex supplemented by boundary 
conditions (45), (48). Provided that the BEP (53) is 
satisfied, the aforementioned complement can 
be obtained in the form of the generalized potential 
(50), uy ei = P ex £ 7 . Difference boundary projection 
(52) and difference BEP (53) are analogous to the 
continuous boundary projection (8) and continuous 
BEP (9), respectively. 

Recall, we have formulated the periodic AP so 
that its solution approaches the solution of the 
infinite-domain AP on any finite neighborhood of T 
as the period Y increases. This gives us grounds (see 
Section 1) for setting the continuous ABC’s in the 
form of BEP (9), in which the Calderon projection 
Pr is calculated using the new finite-domain rather 
than the original infinite-domain AP. Here, we, in 
turn, approximate the solution of the continuous pe- 
riodic AP by the solution of the finite- difference pe- 
riodic AP, which presents the next key step of the 
entire procedure (we mean the next one after in- 
troducing the periodic formulation). This two-step 
scheme leads us to a somewhat non-standard con- 
cept of convergence for the solutions of the differ- 
ence AP. Namely, we will consider convergence of 
the difference solution to the solution of the infinite- 
domain continuous AP on some finite fixed domain 
(e.g., on any rectangle \x — a\ < a, |y| < a, where 
a < X/2, a < Y f 2) as (h x , h y , Y) — - (0, 0, +oo). 
We have already discussed the reasons and conse- 
quences of considering the convergence on a fixed fi- 
nite subdomain only. We should also emphasize that 
the convergence is considered not only as the grid 
size vanishes but also as the period Y synchronously 
grows. Note, to achieve some initially prescribed ac- 
curacy, one should increase the period and decrease 
the grid size consistently. Some estimates connecting 
the grid size, the period, and the desired accuracy 
can be found in ^ . 

Following the considerations of Section 1, one 
can conclude that the convergence of the forego- 
ing type is sufficient for the purpose of construct- 
ing the ABC’s. We will, therefore, use the difference 


analogues (50), (52), and (53) to the continues po- 
tentials (6), projections (8), and BEP’s (9), respec- 
tively, to set the ABC’s in the discrete framework. 
Later on, we will comment on how to choose the 
specific values of the grid size and the period; in 
practice, this choice is always done on the basis of 
the numerical experience. 


2.4 An Application of the Difference 
Potentials and BEP’s for Setting 
the ABC’s 

Let us denote by v the set of those nodes of the 
C-grid that actually determine the penultimate co- 
ordinate line T (see Figure 1); analogously, nodes u\ 
will correspond to Ti. Additionally, we introduce 
on r the set of collocation points to, which is also 
called the collocation grid. This collocation grid will 
be used for specifying the unknowns; typically, it is 
coarser than the grid v. The size of the collocation 
grid u) is not arbitrary, it is connected to the size of 
the Cartesian grid A f°, some relevant estimates can 
be found in ^ . For the practical purposes, we often 
take the total number |w| of nodes of the colloca- 
tion grid u> proportional to the square root of the 
total number |^y | of nodes that constitute the grid 
boundary 7. We should also note that the colloca- 
tion grid is usually not uniform; as a rule, it is more 
concentrated towards the wake region. 

We will approximate the space of clear traces 
£r (see Section 1) by the finite-dimensional space 
Specifically, are the eight-component, vector- 
functions defined at the nodes w, the components of 
contain the trace of the solution u and the trace 
<9u f 3u 

Cw = u, 


of its normal derivative 


We 


V“’ sc 

assume that there is an operation Rrw of interpo- 
lation along T so that as the collocation grid u> is 
refined the continuous function Rrw£w approaches 
the corresponding in the sense of a sufficiently 
strong norm. 

We also introduce the operation of continuation 
of the boundary data from the continuous artifi- 
cial boundary T to the grid boundary 7. Assuming 
that the size of the Cartesian grid Af° is reasonably 
small, one can say that all the nodes 7 are located 
in some small neighborhood of T. Therefore, con- 
sidering £r as the given data, we can drop normal 
from each node 7 to T and then use the first two 
terms of the Taylor expansion to calculate £ 7 . We 
will designate this operation of continuation by 7r 7 r, 
7r 7 r£r = £7. Combining 7r 7 r with the previously in- 
troduced interpolation Rr w , we obtain the operation 
of continuation of the discrete data from T to 7, 
'ft-'/uAuj — 7r 7 rRrw£w — £7 • 
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Finally, recall that the actual data on the curve 
r is calculated numerically at the nodes v. There- 
fore, we will further need the operation R w „ of one- 
dimensional interpolation along T that for a specified 

^ would give £ w , — R wl/ £„. Since 


6 


U; ac 


both and are vector-functions of the same di- 
mension, the interpolation R Wi , is implemented com- 
ponentwise. 

The finite-difference boundary projection P 7 of 
(52) and BEP (53) can be used for setting the ABC’s 
differently. We will begin with the brief descrip- 
tion of the algorithm of 17 , which chronologically 
has come first. 

Substituting the expression £ 7 — 7r yuJ ^ w into the 


BEP (53), we obtain the following equation 
(I 7 P7) ’K'yujfiu! — 0 w 


(54) 


(see 37 ) and therefore 


with respect to Equation (54) can be treated as a 
certain implicit relation between the solution u and 

du 

its normal derivative — at the nodes u>. For the dif- 

?C . 

ference boundary projection P 7 obtained on the ba- 
sis of the central-difference approximation of system 
(40), we actually solve equation (54) with respect to 

1 . . . <9u 

the normal derivatives - 7 — 

.0C. 

express the normal derivatives explicitly in terms of 
u w . The method we employ for solving equation (54) 
is based on the application of a certain variational 
approach, see for more details. Note that in the 
literature, the operators that express boundary val- 
ues of normal derivatives of the solution to a PDE or 
system of PDE’s in terms of boundary values of the 
solution itself are called the Poincare-Steklov oper- 
ators 33 or Dirichlet-to-Neumann maps 2> 3 . 

Having obtained the discrete Dirichlet- 
to-Neumann map S w by solving (54), we are then 
able to calculate £ 7 for any Uj, provided from inside 
the computational domain D{ n : 


£7 — 7*7 w (Rwi/Uy, S w R alJ/ ll^) . (5b) 

Finally, we use the function £ 7 of (55) as the den- 
sity of the generalized difference potential (50) and 
interpolate this potential from the grid J\f ex to the 
nodes ^ C Tj using some local formulas of suffi- 
ciently high order 17 . In so doing, we obtain the 
desirable ABC’s in the form 


U-j/j — R// lJ V’ ex Pe^TT'yw (Rojj/Uj,, S^iR^i/U^) — (56) 

T u„, 

where R^jy/;* is the aforementioned interpolation 
from J\f & x Boundary conditions (56) obviously 


close the discrete system solved inside Di n because 
they provide for the missing relations between the 
values of the solution at the penultimate and out- 
ermost coordinate rows of the C-grid. Moreover, 
we are guaranteed that this closure is right from 
the standpoint of the far-field asymptotic behavior 
of the solution since the operator T of (56) is con- 
structed using the resolved form S w of the boundary 
projection and the potential (50). We should also 
note that since all the operators involved in (56) 
are linear we can actually calculate the matrix of 
the operator T in some appropriately chosen basis. 
This makes the practical implementation of bound- 
ary conditions (56) particularly easy even in spite of 
their nonlocal nature because this implementation 
is, in fact, reduced to a matrix-vector multiplication. 
Note, the simplicity of practical implementation of 
the DPM-based ABC’s is not affected by the shape 
of artificial boundary T ; the matrix form (56) of the 
ABC’s always remains the same although for the 
different shapes of T and F i the matrices T are also 
different. Moreover, as could be seen from our pre- 
vious considerations these different matrices T are 
themselves calculated by means of one and the same 
(i.e., geometrically universal) numerical algorithm, 
which simply uses the shape of the artificial bound- 
ary (more precisely, the actual locations of nodes u 
and u i) as the input data. This algorithm requires 
one solution of the difference AP per basis vector, 
as well as some special numerical procedure that in- 
cludes matrix inversion and multiplication for ob- 
taining S w . Note, the basis, in which we actually 
calculate the matrices of all operators, is actually 
chosen in the space of £ w ’s, and the interpolation 
R wi , is applied afterwards. It is also important to 
emphasize that the RHS’s for the difference AP are 
always concentrated near T (see (51)), and the so- 
lution of the AP also needs to be known only near 
T and Ti. Therefore, the solution of the difference 
AP (direct and inverse Fourier transforms and the 
solution of systems (46), (48) for all k ) requires only 
0(M • J) floating-point operations; a detailed justifi- 
cation of this estimate is contained in 17 , see also 3( ^. 

Further delineation of the foregoing numerical al- 
gorithm for calculating the matrix T from (56) can 
be found in our work . The results of implementa- 
tion of boundary conditions (56) for flow computa- 
tions are reported in 25, 26 , some of these results are 
reproduced and discussed in Section 3 of this paper. 

Another approach to setting the DPM-based 
ABC’s is based on the direct implementation of 
boundary projections, this approach has recently 

r\rr 70 -1 /? ^ 

been proposed in , see also 10 . The main pur- 

pose of introducing the new approach was to reduce 
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the computational cost of the ABC’s. In the new 
methodology, the difference boundary projection P 7 
of (52) is constructed on the basis of the thin-layer 
system (41) using the operator L defined in (43); 
the operation 7r 7w that we use for the continuation 
of boundary data from T to 7 is the same as de- 
scribed above (it is based on the Taylor expansion). 
However, unlike the previous case 17, 25, 26 , in which 
the boundary conditions are driven only by (see 

(9u 

(56)), we now consider both u„ and — as the in- 

oQ h 

put data for the ABC’s. Assuming that these data 
are provided from inside the computational domain 
Di n , we then obtain the density of the generalized 
difference potential by first continuing the boundary 
data from Y to 7 and then applying the difference 
projection P 7 , 


£7 — P 7 TT'yu; P-ui 1 / | U 


du 

’ 5C 


(57) 


In other words, we project the arbitrary boundary 
data provided from inside D in onto the “right man- 
ifold” in the sense that £ 7 of (57) already belongs 
to the image of the boundary projection P 7 , £ 7 E 
ImP 7 , and can therefore be represented as a trace of 
some Ujv^ that solves the equation EhUj^f ex = 0^ er 
and satisfies boundary conditions (45), (48). 

To actually find the aforementioned complement 
uaA*> of£ 7 on J\f ex , we need to compute the difference 
potential P ex £ 7 (see (50)) for the density £ 7 of (57) 
(which requires the solution of the difference AP). 
Then, interpolating this potential from J\f ex to U\ C 
r 1 , we obtain . Combining the foregoing steps, 
we can write the desirable ABC’s in the form 


u 




^iAC e iP ea? P7 TT-yu’ R-o; v 



(58) 


Let us now consider an arbitrary function £ 7 . By 
definition (see formula (50)), the generalized poten- 
tial P eiP (; 7 with the density £ 7 satisfies on J\f ex the 
equation EhPex^-y = ®M ex and boundary conditions 
of the AP (45), (48). In turn, one can easily show 
(see 7 ) that any function ua/;, c that solves the equa- 
tion LftUjy e;c = 0 M ex with boundary conditions (45), 
(48) can be represented as u .j^ ex = P ex r}y, where 
i) y — Tr 7 u J v' e;c . In our case r; 7 — Tr 7 P ex .£ 7 = P 7 £ 7 
and theiefore Pecc^' — ^LV e x — P exV'y — 

In other words, for any £ 7 the following relation 


P ex^-y — Pea;P 7^7 


is true. Consequently, instead of (58) we can write 
the ABC’s as follows 


e.x'K-y 


U, 


T 




(59) 


(the matrices T in equations (56) and (59) are obvi- 
ously different). Clearly, boundary conditions (59) 
provide for the missing relations between the values 
of the solution on the penultimate and outermost 
coordinate rows of the C-grid and therefore close 
the finite-difference system that we solve inside D in . 
Moreover, this closure is consistent with the desir- 
able far-field behavior of the solution because of the 
projection P 7 incorporated in (58). 

Comparing boundary conditions (56) and (59) we 
see that they are essentially different. Direct usage 
of the boundary projection P 7 in (58), (59) allows us 
to completely avoid the entire resolving stage of the 
algorithm that is inherent for the approach of the 
first, type summarized in formula (56). Recall, the 
resolving stage in (56) is associated with the compu- 
tation of (resolved form of the boundary projec- 
tion) on the basis of equation (54). Elimination of 
this stage implies an essential simplification of the al- 
gorithm, as well as noticeable economy of computer 
resources, i.e., the reduction of the computational 
cost of boundary conditions (59) in comparison with 
boundary conditions (56). 

Of course, another part of this cost reduction, 
which is even more essential, is accounted for by 
the reduction of order n of one-dimensional finite- 
difference system (46) from the eighth to the fourth. 
Indeed, we recall that boundary conditions (56) were 
obtained on the basis of the second-order central- 
difference discretization of system (40) and bound- 
ary conditions (59) were obtained on the basis of 
the discretization (43) of system (41). In so do- 
ing, the matrices A& and in system (46) have 
order n — 8 for boundary conditions (56) (see ^ 7 ) 
and order n — 4 for boundary conditions (59) (see 
(44), (47)). According to the solution of one- 
dimensional problem (46), (48) for each k costs 
G(M • n 2 ) floating-point operations. Therefore, the 
solution of the entire AP costs 0{M ■ J ■ n 2 ) opera- 
tions. For both boundary conditions (56) and (59) 
we have to solve the AP repeatedly, one time per 
basis vector. Consequently, one can expect that for 
the same geometry of the discrete sets v, , and uj, 
for the same basis in the space of £ w ’s, and for the 
same grid Af° , the computational cost of the ma- 
trix T from (59) will be at least four times less that 
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the cost of matrix T from (56). Taking into account 
the foregoing elimination of S w from the structure of 
ABC’s, we conclude that the overall improvement of 
the computational efficacy when going from (56) to 
(59) will be even more drastic. Our computational 
experiments do corroborate this theoretical expec- 
tations. In fact, we could gain up to a factor of 
five in the reduction of cost of the ABC’s (56) in 
comparison with ABC’s (59). Moreover, our prelim- 
inary estimates show that in the case of three space 
dimensions this gain may increase. 

As mentioned above, along with being computa- 
tionally cheaper than the original technique (56) the 
new methodology (59) also appears much simpler 
from the algorithmic standpoint. At the same time, 
boundary conditions (59) do posses all the afore- 
mentioned favorable properties that are relevant to 
boundary conditions (56). Namely, they are geomet- 
rically universal, easy to implement in practice, and 
of course, they perform as good as (or even better 
than) (56) from the standpoints of accuracy, overall 
efficacy, and robustness (see Section 3). 

2.5 Implementation of DPM-based ABC’s 

The ABC’s of both types (56) and (59) are de- 
signed to close the finite-difference system solved 
inside Di n , i.e., to make sure that the number of 
equations solved inside Di n is equal to the number 
of unknowns. Both relations (56) and (59) are spa- 
tially nonlocal, which in practical terms means that 
the matrices T are dense. Although these matri- 
ces are, in fact, structural, we have not used this 
property for practical computations yet. (Qualita- 
tive structure of the matrices T can be understood 
from physical considerations. If the vectors u„ and 
u Vl are arranged properly, then we can consider T 
as being composed of several blocks. Each block of 
T would correspond to one physical variable (u, v, 
p , or p ) row-wise, i.e., for all nodes v, and one phys- 
ical variable column-wise, i.e., for all nodes V\, and 
would have a kind of “diagonal dominance” in the 
sense that the entries located near the main diago- 
nal will be greater that those located far away from 
this diagonal. In physical terms, it merely means 
that each specific node influences its close neighbors 
stronger than it influences the nodes located on the 
other side of the computational domain.) 

So far, we have been discussing the ABC’s only 
from the viewpoint of closing the system solved in- 
side Di n so that the closure is consistent with the 
desirable far-field behavior of the solution. In prac- 
tice, however, the construction of a formal closure of 
the finite-difference system solved inside Z>j n is not 
sufficient, we also have to combine this closing proce- 


dure with the specific solver. The majority of solvers 
currently used in CFD for calculating the steady- 
state viscous flows on the basis of finite-difference 
discretizations employ various types of pseudo-time 
iterations. In most cases, the iterations are explicit 
in time and may be enhanced by different techniques 
for the purpose of accelerating the convergence, e.g., 
by multigrid. 

We have to emphasize that both boundary condi- 
tions (56) and (59) are particularly well fitted for the 
combined usage with explicit iterative solvers. In- 
deed, let us assume that on some time level (which, 
in particular, may be zero) the solution is already 
known on the entire C-grid. Then, advancing one 
time step by means of some explicit technique we 
obviously cannot obtain the next-level solution also 
on the entire C-grid. The reason for that is exactly 
the same as why the original steady-state finite- 
difference system inside D{ n would be subdefinite 
without the ABC’s — the stencil applied to some 
external nodes of the C-grid may partially “fall out” 
of the domain. In other words, when using solely 
the procedure employed inside Di n , we can obtain 
the solution on the upper time level only at the “in- 
ternal” nodes of the C-grid. For the case of the 
stencil 3x3, this “internal” set includes all nodes 
of the C-grid except for the outermost coordinate 
row v i C Ti. The values of the solution on this out- 
ermost coordinate row should therefore be provided 
by the ABC’s so that the solution on the upper time 
level becomes available everywhere, which makes the 
next iteration feasible. Looking at boundary condi- 
tions (56) and (59) one can see that the desirable 
complement of the solution on the upper time level 
can easily be obtained by means of either one of 
these techniques. When using boundary conditions 
(56), we simply take that is already computed 
on the upper time level and, applying the operator 
T (matrix- vector multiplication), obtain u^. This 
operation is repeated on every iteration; if the re- 
laxation procedure requires evaluation of the resid- 
uals more than once per iteration (e.g., multi-stage 
Runge-Kutta), then boundary conditions (56) are 
used as many times per iteration as the residuals 
need to be evaluated. 


The implementation of boundary conditions (59) 
is analogous. After making one iteration of the 
Navier-Stokes solver inside Di n we know the solu- 
tion on the upper time level everywhere in the inte- 
rior of the curve T. Therefore, we can take on 
the upper time level as done above and also can eas- 
du | 

using the available data. Then, 


ily calculate 


dC 


applying the matrix T from (59), we obtain u (/1 on 
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the upper time level and therefore make it possible to 
advance another time step. One can see that in both 
cases (56) and (59) practical implementation of the 
DPM-based ABC’s is very easy since it is reduced 
to a matrix- vector multiplication on each iteration. 
Moreover, the implementation is in no way affected 
by the shape of artificial boundary, although the op- 
erators T themselves would, of course, depend on 
the geometry. 

The implementation of the DPM-based boundary 
conditions would obviously require some changes if 
the solver employed inside Di n uses multigrid for the 
acceleration of convergence. Namely, for the multi- 
grid iterative solver the values of the solution on the 
outermost coordinate row of the grid should be pro- 
vided every time the residuals need to be evaluated 
on every level of multigrid. Since the operators T 
do depend on the geometry, we may formally need 
to calculate a separate operator for each subsequent 
grid. It, however, turns out that at least for some 
multigrid strategies the numerical process appears 
to be sensitive only to the ABC’s specified at the 
finest level of multigrid, whereas the sensitivity of 
the numerical process to the boundary conditions 
on all the coarser levels is negligible. 

In all the computations reported in Section 3, we 
have used the algorithm 33 , 34 g wanson anc j 
Turkel for obtaining the steady-state solutions of the 
Navier-Stokes equations inside D jn . This algorithm 
is based on the second-order central-difference ap- 
proximation in space and Runge-Kutta relaxation 
in time. In practice, we have always used five-stage 
Runge-Kutta time stepping. The code, in which the 
algorithm of ’ ’ is implemented, allows one to 

use different multigrid strategies for the acceleration 
of convergence. Depending on the specific computa- 
tional variant, we used from three to five levels of 
multigrid with W-cycles, when one iteration is done 
on the finest level and two iterations are done on 
each of the coarser levels of multigrid. As the addi- 
tional means of convergence acceleration, the algo- 
rithm of ’ ■ includes local time stepping and 

residual smoothing. 

The standard treatment of the external boundary 
in the code ’ ■ is based on the locally-one- 

dimensional analysis of characteristics for the inflow 
part of the boundary and the boundary conditions 
of extrapolation for the outflow part of the bound- 
ary. This approach is actually one of the most well- 
known and widely used in CFD, its different ver- 
sions have many times been described in the litera- 
ture, see, e.g., the reviews ’ ’ . The purely local 

characteristics/extrapolation treatment may or may 

qo 09 9 A 

not be enhanced in the code ’ ’ by the point- 


vortex correction. This lift-based correction ^ usu- 
ally improves the results provided by the original 
local ABC’s. 

We have experimentally made certain that for 
the different flow regimes (see Section 3, as well 
as ^ 6 ) computed on the basis of the foregoing 
multigrid strategy, neither the convergence history 
nor the calculated solution depend on whether we 
specify the standard local ABC’s (see above) on all 
levels of multigrid or we specify these boundary con- 
ditions on the finest level only and for all the coarser 
levels simply retain the boundary values provided 
from the finest level. This gave us reasons to ex- 
pect that the nonlocal DPM-based ABC’s (56) or 
(59) can also be set on the finest level of multigrid 
only, whereas the boundary values for all the coarser 
levels will be provided from the finest one. In all 
the computations reported below, we used exactly 
this scheme of implementation of the nonlocal DPM- 
based ABC’s. Numerical results (see Section 3) cor- 
roborate the possibility of doing so, at least for those 
multigrid strategies that we used for our computa- 
tions. 

Finally, we should note that the implementation of 
the DPM-based ABC’s with implicit iterative solvers 
also seems feasible, at least from the formal stand- 
point. Although we have never tried to run any nu- 
merical experiments, theoretically this implementa- 
tion would simply mean that the nonlocal relations 
of type (56) or (59) are incorporated in the system 
that is solved on the upper time level on each step 
of the implicit iteration process. 

2.6 Miscellaneous Issues Important 

for Calculation of DPM-based ABC’s 

First, we comment here on the choice of the dis- 
cretization parameters for the difference AP. The un- 
knowns for this problem are specified on the Carte- 
sian grid A" 0 . The cell size of this grid should be 
chosen so that the distance between the curves F 
and Ti (see Figure 1) is resolved. In practical com- 
putations we usually take the cell size of the grid J\f° 
to be 3-5 times smaller than the minimal distance 
between T and Fi. 

As concerns another important parameter of the 
difference AP, the period Y (see Figure 1), it should 
be chosen so that to ensure the sufficient accuracy 
of computations. This choice, of course, cannot be 
done without taking into account the previous com- 
putational experience. However, we first have to 
choose the unit for measuring Y. It is reasonable 
to expect that the average diameter of the compu- 
tational domain D{ n would make a better unit for 
measuring the period Y than the characteristic size 
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of the immersed body(ies) (e.g., airfoil chord, see 
Figure 1). Indeed, the final accuracy will depend on 
the “extent of convergence”, i.e., on how close the 
solution of the periodic AP has approached the solu- 
tion of the original nonperiodic AP. Since the size of 
the computational domain Di n actually determines 
the size of that subdomain in Dy , on which we con- 
sider the convergence (e.g., the rectangle |y| < o, 
\x — a\ < a), then exactly the size of Di n becomes 
a natural unit for measuring the period Y as the 
variable that controls the convergence. 

As mentioned above, the value of the period Y 
is one of the factors that determine the accuracy of 
computations. The accuracy, of course, also depends 
on the actual size of Di n . Generally, the DPM- 
based ABC’s allow one to use much smaller compu- 
tational domains than other available methods do. 
The specific numerical results will be reported in 
Section 3; here, we provide only for a qualitative pic- 
ture. In ^ , we have conducted a series of compu- 
tations for the relatively small values of the period 
Y , about 2-4 diameters of Di n . It turns out that for 
a small computational domain (2-3 chords of the 
airfoil for transonic flow regions), the solution ob- 
tained on the basis of the DPM-based ABC’s differs 
within 2-4% from the asymptotic solution obtained 
in a very large computational domain (about 30-50 
chords depending on the specific variant). (Note, 
when comparing the two solutions, we actually com- 
pare the corresponding force coefficients: lift, drag, 
skin friction). These results (see > 2 ®) can be satis- 
factory for some cases; however, when the accuracy 
needs to be improved one has to choose larger pe- 
riods Y. For example, to keep the force coefficients 
within 0.01%— 0. 1% of the corresponding asymptotic 
values, one may need to choose the period Y of about 
50 diameters of the computational domain or more. 

That big periods Y , along with the small cell size 
of the grid J\f° prescribed by the distance between F 
and Ti may, generally speaking, result in an expen- 
sive numerical procedure for computation of the op- 
erators T. To reduce this cost, we introduce nonuni- 
form grids with respect to y. For example, the cell 
size of such a grid may remain constant (the same 
as it would be for the uniform grid) in the neigh- 
borhood of Di n and then enlarge as |y| increases. 
Clearly, the discrete Fourier transform in the y di- 
rection that we have formerly used for the separation 
of variables does not apply to the stretched grids. 
Consequently, we will need some other procedure to 
separate the variables and to therefore reduce the 
difference AP to a family of one- dimensional sys- 
tems (46) with boundary conditions (48). The sepa- 
ration of variables for a finite-difference counterpart 


of system (41) means that the discrete transform 
should simultaneously diagonalize the difference ap- 
proximations to both the first and the second deriva- 
tives with respect to y. We have not studied thor- 
oughly the question of whether or not one can find 
a special distribution of nodes in the y direction and 
some consistent with this distribution discretization 
so that the diagonalizing transform appears orthog- 
onal. In practice, it turns out that usage of the skew 
bases solves the problem of simultaneous diagonal- 
ization of the first and the second derivatives and 
at the same time allows one to still maintain high 
accuracy of the final results. 

Following this idea we first introduce some second 
order discretization of the first derivative with re- 
spect to y on the stretched grid, we also take into 
account the periodicity conditions (45). In the ma- 
trix form, this operation can be represented as a cer- 
tain (2 J + 1) x (2 J + 1) matrix D, which has three 
non-zero diagonals, j — 1 < i < j + 1, and also non- 
zero off-diagonal corner entries. Then, the matrix 
V 2 represents the second difference derivative with 
respect to y. Both these matrices can obviously be 
diagonalized by the same transform, which is com- 
puted in practice with the help of the standard li- 
brary eigenvalues/eigenvectors subroutines (IMSL). 
The rest of the algorithm remains the same as de- 
scribed above. The only difference is that instead of 

and t & in (47) one should plug in the eigenvalues of 
V and their squares, respectively. Our experiments 
show that usage of the stretched grids can drastically 
reduce the computer effort required for calculating 
the nonlocal DPM-based ABC’s. An essential part 
of the results reported in Section 3 is obtained on 
the basis of the stretched-grid algorithm. 

Another means for reducing the computational 
cost of the DPM-based ABC’s is implementation of 
the algorithm on parallel platforms. There are, gen- 
erally speaking, a few ways to parallelize the compu- 
tation of operators T. Since we anyway calculate the 
matrix of the linear operator in a certain basis, it is 
possible to calculate independently and, therefore, 
simultaneously, the columns of this matrix, where 
each column corresponds to its own basis vector. 
We, however, have chosen another way of paralleliza- 
tion, which, in our opinion, requires minimal mod- 
ifications of the already developed sequential code. 
Namely, one-dimensional systems (46) with bound- 
ary conditions (48) are obviously independent be- 
cause they are obtained by the separation of vari- 
ables ( k is a parameter). Therefore, these systems 
can be solved in parallel on the different processors 
of a multi-processor computer. This approach has 
been implemented on an eight-processor CRAY Y- 



DIFFERENCE POTENTIALS METHOD FOR EXTERNAL CFD 


21 


MP. Numerical experiments show the reduction of 
the “wallclock” time required for calculating the op- 
erator T of (59) by up to a factor of five in compari- 
son with the standard single-processor implementa- 
tion. 

Finally, we should note that the treatment of tur- 
bulence in the far field requires special attention. In- 
deed, the parameter Re in systems (40) and (41) rep- 
resents the true molecular Reynolds number, there- 
fore these forms of the governing equations apply 
directly only to the laminar flows. The approximate 
treatment of turbulence in the far field for the pur- 
pose of constructing the ABC’s was proposed in the 
work ^ . Below, we reproduce some results of this 
work. 

To numerically simulate turbulent flows, we use 
an algebraic turbulence model (Baldwin-Lomax) in- 
corporated in the code 3 ^’ 33, 3 ^. This model is rel- 
evant to describing the flow in the vicinity of the 
immersed body(ies). In the far field, we use simpler 
approach based on the concept of effective turbulent 
viscosity 3 ®. The idea is to qualitatively describe 
turbulent flow (i.e., the process of turbulent mixing) 
as a laminar flow of model fluid having some new 
“turbulent” viscosity. 

To obtain the relation between the molecular and 
effective turbulent viscosity, we use the following 
considerations. First, we refer to the incompress- 
ible case and consider laminar flow. Here, we have 
the following distribution of u-velocity (perturba- 
tion with respect to the far-field value t<o) 36 : 


W 

u = — = exp 
2\Ar p 2 vupx 



(60) 


We recall 36, 37 that formula (60) is obtained under 
the natural assumption that the far-field solution ac- 
tually depends neither on the shape of the immersed 
body nor on the type of flow in its close vicinity but 
only on one constant W, which is the total drag. 
(Note that the dimension of W in (60) is that of 
force per unit length.) 


For the turbulent case we assume 36, 37 that the 
mixing length / for the wake flow is proportional to 
the local width b of the wake, l/b = const. Then, ap- 
proximately replacing the value of the derivative in 

i „ . _ . . l0 \du 

the expression for turbulent viscosity, v t = / 


dy 

by the ratio u max /b, where u max is the maximal de- 
viation of the actual velocity from uq (which corre- 
sponds to the middle of the wake), we easily obtain 
u t — const b -Umax, which, in particular, implies that 

Of? 

v t does not depend on y. We may also assume 
that analogously to the laminar case the wake-type 


solution for the turbulent flow is self-similar. Then, 
we have u = u max -(j)(y/b), which immediately yields, 


W 


= puo i u dy — pupbu max J <j> (j) d (jf) 
-b -1 


and consequently, 


kW 


Ut PUQ ’ 


(61) 


where k = const. Once can easily see from (61) that 
u t = const throughout the whole wake region and 
therefore, the structure of turbulent wake behind the 
body (i.e., the profile of mean velocity) appears to 
be the same as the structure of laminar wake 3 ®’ 37 
since we actually have the same solution as (60), 


W 


u = 


VtU pX 

1 fw 


exp 



exp — 


u 0 y 2 \ _ 
4 u t x ) 

puo V\ 

4 kWx ) ’ 


(62) 


only for some other constant u t instead of the true 
molecular viscosity v. 

Recall, our purpose is to find such u t that, being 
substituted into (60), will provide the same wake 
solution as we would have for the turbulent case. 
Clearly, v t from (61) satisfies this requirement (see 
(62)), so, let us now determine the specific value of 
k. Since the effective viscosity is constant through- 
out the whole wake region, we assume that in the 
far field it preserves the same value as it has in the 
outer part of the boundary layer near the trailing 
edge of the immersed airfoil. Using the Clauser con- 
jecture 3 ® to calculate the latter quantity, and re- 
stricting ourselves (for qualitative consideration) by 
the case of not high longitudinal pressure gradients, 
we can obtain, 


k c W 

vt = , 

pu 0 

which means that the value of the unknown constant 
k in (61) may be chosen the same as the value of the 
Clauser constant, kc — 0.0168. 

To independently determine W, we recall that in 
the case of turbulent flow past a flat plate the drag 

Of? ^ 

is given by 

W = O.QmRe- 1 / 7 pu 2 0 L, (63) 


where Re is the actual Reynolds number based on 
the molecular viscosity and L is the characteristic 
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length. To take into account the compressibility, we 


multiply equation (63) by 


2 + 0.5(7 “ 1) Mq 


5/7 


in accordance with the Tucker conjecture 36 . Then, 
we introduce an empirical constant 0, which is the 
ratio of the pressure and viscous contributions to 
the total drag (obviously, 0 = 0 for the flat plate); 
for the transonic flows computed below, 0 & 5/2. 
Finally, we obtain the following relation, 


Re t = 


uqL 

Vt 


Re\ /7 ( 2 

O.O3O7fc c (0 + 1) \2 + 0.5(7 ~ 1) Mq 

which is used for determining the effective turbulent 
Reynolds number in all the turbulent computations 
presented in the next section. We emphasize that 
the proposed treatment of turbulence in the far field 
is only qualitative. However, this approach seems to 
be be justified by the numerical experiments. 



3 Numerical Results 


3.1 Acceleration of Convergence 

One of the most important aspects of implementa- 
tion of any ABC’s is the influence that the boundary 
conditions exert on the convergence to steady state. 
Our numerical experiments for both two 2 ^' 26 anc j 
three 38, 39 space dimensions show that the nonlocal 
DPM-based ABC’s can essentially speed up the con- 
vergence of the multigrid iterations compared to the 
standard characteristics-based boundary conditions. 

This positive influence on the convergence rate is, 
however, not a general situation. It appears that 
the acceleration of convergence typically occurs only 
when the interior iterative solver involves multigrid. 
Otherwise, the nonlocal highly accurate ABC’s ei- 
ther do not influence the convergence at all or may 
even slow it down. For example, the observation of 
the latter kind was done by Ferm in 38 for the non- 
local ABC’s 38, 39 (by Gustafsson and Ferm) that 
are constructed for the Euler flows in ducts using 
Fourier transform in the cross-stream direction. The 
same phenomenon also occurs for the external in- 
viscid flows as shown by Ferm in the work 49 , in 
which he studies the convergence of pseudo-time it- 
erations for the Euler equations supplemented by 
the nonlocal ABC’s 43 (boundary conditions 43 are 
constructed analogously to 38 ’ 39 for elliptic artifi- 
cial boundaries). To accelerate the convergence of 
pseudo-time iterations with nonlocal boundary con- 
ditions, Ferm in 38 ’ 49 employs the technique of 42 
by Engquist and Halpern (or its modification), which 


allows him to make the convergence at least as fast 
as it is for the simplest locally-one-dimensional non- 
reflecting boundary conditions that are based on the 
analysis of characteristics. On the other hand, when 
boundary conditions 43 are implemented along with 
some multigrid Euler solver inside the computational 
domain they no longer slow down the convergence 
and, therefore, no longer require special acceleration 
procedures (like the one from 42 ). This has been 
demonstrated by Ferm in the work 43 , in which he 
shows that in order to reduce the initial error by a 
prescribed factor one needs roughly the same num- 
ber of multigrid cycles for both nonlocal ABC’s 43 
and the characteristic non-reflecting boundary con- 
ditions. 

As mentioned above, when implemented along 
with the multigrid algorithm 32 ' 33 1 34 (see Sec- 
tion 2), the DPM-based ABC’s are capable of even 
speeding up the convergence to steady state com- 
pared to the standard boundary conditions. Let 
us reproduce here several graphs from 2 ^ that rep- 
resent the convergence history for different sub- 
sonic and transonic laminar flows around the airfoil 
NACA0012. In the captions to the figures below, a 
denotes the angle of attack. 

From Figures 2, 3, 4, and 5, one can easily see 
that usage of the DPM-based ABC’s can increase 
the convergence rate of the multigrid iterations by 
up to a factor of three depending on the specific 
variant of computations. Note, the subcritical (i.e., 
fully subsonic) laminar cases that correspond to Fig- 
ures 2, 3, 4, and 5 have been computed on the grids 
with small stretching ratios because near the air- 
foil surface those grids could be chosen relatively 
coarse. As a result, we have used global rather 
than local Courant step for iterations in time. In 
this respect, one can say that Figures 2, 3, 4, and 
5 demonstrate the influence exerted by the DPM- 
based ABC’s on a “pure” multigrid (augmented only 
by residual smoothing). 

For the case of two-dimensional turbulent flows 
that are computed on the grids with much higher 
stretching ratio and with the local (i.e., chosen cell 
by cell) Courant step in time, we have not been 
able to obtain as drastic convergence speedup as 
for the foregoing laminar cases. The history of con- 
vergence for two different two-dimensional transonic 
turbulent cases is presented in Figures 6 and 7. We 
however, mention, that for many three-dimensional 
transonic turbulent cases, the DPM-based ABC’s 
have been able to produce the increase of the con- 
vergence rate about as big as shown above for the 
two-dimensional laminar flows. The corresponding 
results are reported in our work 38, 19 and will also 
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Figure 6. Convergence history: log|jp res ^ U a?||oo versus number of cycles; RAE2822, Mo = 0.73, a — 2.79°, 
Re = 6.5 • 10 6 . Average radius of Di n about 11 chords, normal spacing near the airfoil 10 -4 . 
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Figure 7. Convergence history: Iog||/VesieZua/||oo versus number of cycles; RAE2822, M 0 = 0.73, a = 2.79°, 
Re = 6.5 • 10 b . Average radius of Dj n about 11 chords, normal spacing near the airfoil 10~ 5 . 
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Figure 8. Convergence history: logH/Ve^ua/lloo versus number of cycles; NACA0012, Mo = 0.85, a = 1°, 
Re = 4000. Grid 256 x 64, average radius of Di n about. 5.5 chords. 
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Figure 9. Convergence history: log||/9 res? rf ua ;jjoo versus number of cycles; RAE2822, M 0 = 0.73, a = 2.79°, 
Rt = 6.5 • 10 6 . Average radius of D{ n about 6 chords, normal spacing near the airfoil .33 • 10 -4 . 
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be briefly discussed later in this paper. 

Returning to Figures 2, 3, 4, and 5, we see 
that the convergence rates for two different types 
of ABC’s are the same on the initial stage of the 
iteration process; then, for the DPM-based ABC’s 
the convergence rate remains the same all the time 
and for the standard boundary conditions it drasti- 
cally decreases. Therefore, it would be reasonable 
to assume that the ABC’s start to actually influ- 
ence the convergence only after the numerical per- 
turbations caused by the immersed body reach the 
external boundary. In other words, the DPM-based 
ABC’s become most effective from the standpoint of 
convergence acceleration on the so-called asymptotic 
stage of the multigrid. The similar type of behavior 
can be observed for the three-dimensional computa- 
tions as well (see below). 

We, however, have to say that although the ac- 
celeration of multigrid convergence provided by the 
DPM-based ABC’s is extremely important for ap- 
plications, the mechanism of interaction of the non- 
local DPM-based ABC’s with multigrid may require 
an additional study. In fact, neither rigorous mathe- 
matical explanation of the convergence speedup nor 
a definite experimental conclusion of why and when 
it happens is available as of yet. For our two- and 
three-dimensional computations, we have used dif- 
ferent multigrid strategies (W and V cycles, respec- 
tively), also all the three-dimensional cases have 
been computed with the local time step, and the 
results are also different. Whereas for the two- 
dimensional transonic turbulent flows we did not see 
much of an increase in the convergence rate, in three 
dimensions the strongest speedup occurs right for 
the transonic turbulent cases (see the work 19 
and also below). At the same time, the convergence 
rates for the subsonic turbulent flows in three space 
dimensions are the same for the ABC’s of different 
types 19 . As for the laminar flows, the experi- 
ments have been conducted in two space dimensions 
only. 

Analyzing the influence that nonlocal boundary 
conditions may exert on the convergence of multi- 
grid iterations, we should also note that many mod- 
ern multigrid solvers are not optimal themselves. A 
massive effort is currently underway towards con- 
structing the new finite-difference schemes, for which 
the convergence characteristics of multigrid meth- 
ods would essentially improve. For example, the 
work in this direction has been done by Ta’asan 44 
and Sidilkover 45 . In 44 , Ta’asan devised an essen- 
tially optimal multigrid solver for the Euler equa- 
tions in subsonic regime. Due to the separate treat- 
ment of the elliptic and advection parts of the sys- 


tem, this approach allows one to achieve in sub- 
sonic regime the convergence rates similar to those 
that can be obtained when solving the full poten- 
tial equation. Sidilkover in 4j proposed the so-called 
genuinely multidimensional high-resolution scheme. 
This scheme has stability properties much superior 
to those that are relevant to the standard meth- 
ods and, therefore, facilitates the construction of a 
very simple and efficient multigrid algorithm (using 
Gauss-Seidel relaxation as a smoother) that would 
apply to the entire range of Mach number. Particu- 
lar efficacy demonstrated by the DPM-based ABC’s 
when implemented in combination with multigrid 
gives us reasons to hope that these boundary condi- 
tions may essentially contribute in developing the 
new generation of effective multigrid-based algo- 
rithms. 

Finally, we should mention that the DPM-based 
ABC’s generally improve the robustness of the en- 
tire numerical procedure. In conducting our com- 
putational experiments we have noticed that 

sometimes the multigrid iterations ^ 4 supple- 

mented by the standard characteristic/extrapolation 
boundary conditions simply fail to converge, which 
never happens if these standard ABC’s are replaced 
by the nonlocal DPM-based boundary conditions. 
In Figures 8 and 9, we show the history of con- 
vergence for the corresponding computations. The 
similar phenomenon has been observed in three di- 
mensions as well. Namely, for a transonic turbulent 
flow with separation (see ^ 9 ), the multigrid iteration 
procedure with standard boundary conditions failed 
to converge, whereas the DPM-based ABC’s have 
still been able to ensure a fast convergence to steady 
state. 

3.2 Accuracy 

Another most important outcome of usage of the 
DPM-based ABC’s is the essential increase of ac- 
curacy that these boundary conditions provide for 
in computing the external viscous flows. Below, we 
compare some numerical results obtained on the ba- 
sis of boundary conditions (59) for a certain tran- 
sonic turbulent flow around the airfoil RAE2822 
with the results obtained for the same flow regime 
on the basis of the standard local ABC’s (character- 
istics/extrapolation) enhanced by the point- vortex 

or 

correction . In Table 1, we present the results 
for three different grids, 640 x 128, 608 x 112, and 
600 x 104 nodes, that correspond to the computa- 
tional domains of the average radii of 50, 8, and 
2.5 chords of the airfoil, respectively. It is impor- 
tant that each subsequent (smaller) grid is obtained 
here by cutting off several external coordinate lines 
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Table 1. Comparison with the point-vortex (p.-v.) model for RAE2822 airfoil; Mq — 0.73; Re o = 6.5 ■ 10 6 ; 
a = 2.79°; normal grid spacing near the airfoil 0.5 • 10~ 5 . 


Domain “radius” 

3 chords 

8 chords 

50 chords 

Grid 

600 x 104 

608 x 112 

640 x 128 

Type of ABC’s 

p.-v. 

(59) 

p.- V. 

(59) 

p.-v. 

(59) 

C t 

0.8653 

0.8591 

0.8624 

0.8589 

0.8603 

0.8593 

relative error 

0.58% 

0.02% 

0.24% 

0.04% 

0% 

0% 

C d x 10 

0.1203 

0.1263 

0.1209 

0.1261 

0.1255 

0.1260 

relative error 

4.14% 

0.24% 

3.67% 

0.08% 

0% 

0% 

C D x 10 

0.1755 

0.1816 

0.1762 

0.1815 

0.1810 

0.1815 

relative error 

3.04% 

0.05% 

2.65% 

0% 

0% 

0% 


Table 2. Comparison with the point-vortex (p.-v.) model for RAE2822 airfoil; Mq = 0.73; Re 0 = 6.5 • 10 6 ; 
a = 2.79°; normal grid spacing near the airfoil 0.5 • 10 -5 . 


Domain “radius” 

2.5 chords 

50 chords 

Grid 

320 x 64 

320 x 64 

640 x 128 

Type of ABC’s 

p.-v. 

(59) 

p.-v. 

(59) 

p.-v. 

(59) 

Ci 

0.8688 

0.8560 

0.8504 

0.8492 

0.8603 

0.8593 

relative error 

2.15% 

0.38% 

1.15% 

1.17% 

0% 

0% 

C d x 10 

0.1123 

0.1259 

0.1260 

0.1265 

0.1255 

0.1260 

relative error 

10.5% 

0.07% 

0.40% 

0.39% 

0% 

0% 

C f x 100 

0.5469 

0.5492 

0.5478 

0.5480 

0.5543 

0.5544 

relative error 

1.34% 

0.94% 

1.17% 

1.15% 

0% 

0% 

C D x 10 

0.1670 

0.1808 

0.1808 

0.1814 

0.1810 

0.1815 

relative error 

7.73% 

0.39% 

0.11% 

0.05% 

0% 

0% 


of the preceding (bigger) grid. This is done in order 
to completely avoid any possible influence that the 
change of the grid near the airfoil surface may exert 
on the solution. 

From Table 1 one can see that the correspond- 
ing asymptotic values of the force coefficients (lift 
Ci, wave drag Cd , total drag Cd), be., the values 
obtained for the large (50 chords) computational 
domain, are very close to one another for the dif- 
ferent types of ABC’s. However, as the artificial 
boundary approaches the airfoil the discrepancy be- 
tween the corresponding values increases, and the 
force coefficients obtained on the basis of boundary 
conditions (59) deviate from their asymptotic values 
much less than the coefficients obtained using local 
ABC’s. In other words, the nonlocal DPM-based 
ABC’s allow one to use much smaller computational 
domains than the standard boundary conditions do 
and to still maintain high accuracy of computations. 
Moreover, from Table 1 one can see that unlike the 
ABC’s (59), which perform well for all coefficients, 


the point-vortex boundary conditions perform much 
better for the lift coefficient Cj than they do for the 
drag coefficients Cd and Cd- This behavior seems 
reasonable since the point-vortex model is a purely 
lift-based treatment and does not take into account 
drag at all. 

In Table 2 we also compare the results obtained 
using the two aforementioned types of ABC’s; how- 
ever, the computations presented in this table were 
conducted on the different grids. One can see 
that boundary conditions (59) outperform the point- 
vortex ABC’s in these cases as well ( Cf in Table 2 
is the skin friction). 

We should also emphasize that the benefit of us- 
ing smaller computational domains and, as a conse- 
quence, smaller grids, i.e., the grids with lesser num- 
ber of nodes, is not only the direct reduction of the 
computational work because of the grid shrinkage 
but also the improvement of convergence because 
the grids may be chosen less stretched. 
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3.3 Entry-Wise Interpolation 

The computational overhead associated with the 
usage of nonlocal ABC’s (56) or (59) consists of two 
parts. The first part is accounted for by the matrix- 
vector multiplications that we do every time we need 
to evaluate the residuals (see Section 2). These op- 
erations add about 1-2% to the total cost of com- 
putations (more precisely, to the cost of the same 
number of multigrid cycles but without the nonlocal 
ABC’s). Generally, we estimate this additional ex- 
pense as low (taking into account the benefits pro- 
vided by the DPM-based ABC’s), and we do not 
think that any special effort towards reducing this 
cost is currently required. Nonetheless, we note that 
there may still be some room for reducing this cost 
by using the multiresolution-based techniques (see 
our work 18, 21 ). Moreover, we should mention that 
the operation of matrix- vector multiplication is, as a 
rule, fully vectorizeable, which provides for an addi- 
tional advantage of using the DPM-based ABC’s on 
the CRAY machines. Furthermore, we expect that 
the new discretizations that are currently under de- 
velopment (see ^) will allow one to use relax- 
ation procedures other than the methods of Runge- 
Kutta type that are most widely used now. As a con- 
sequence, the overhead associated with the matrix- 
vector multiplications (56) or (59) may be reduced 
since, for example, the methods of Gauss-Seidel type 
would require only one such operation per iteration, 
whereas the multi-stage Runge-Kutta methods re- 
quire as many multiplications (56) or (59) as the 
number of stages is. Finally, we mention that the 
DPM-based ABC’s can be implemented directly, i.e., 
the operation T can be computed explicitly (with- 
out first calculating the matrix in a basis) every time 
u| needs to be updated. This strategy has, in fact, 
been chosen for all our three-dimensional computa- 
tions 20, 21 ; the corresponding results will be 

briefly commented on later in this paper. 

If the matrix-based strategy is used, then the sec- 
ond part of the total overhead due to the ABC’s is 
accounted for by the computational cost of the oper- 
ators T themselves (see (56) and (59)). In the case 
of two space dimensions, we could always keep the 
corresponding cost at a level of about 10% of the 
total work required to calculate the steady-state so- 
lution using the algorithm 32, 33 ’ 34 with the specific 
multigrid strategy discussed in Section 2. Basically, 
this additional expense can be regarded low as well. 

As mentioned above, in three space dimensions we 
have so far been using another strategy that did not 
require the calculation of matrices T at all. How- 
ever, in the future we may need it, especially in 
the view of possible multiple runs. Our prelimi- 


nary estimates show that in the case of three space 
dimensions the cost of the operator T may appear 
higher (in relative terms) than it is for two space di- 
mensions. Therefore, we propose a special approach 
to decreasing the computational cost of the nonlo- 
cal DPM-based ABC’s in the framework of massive 
computations. 

Very often in CFD, one needs to calculate differ- 
ent flows around the same configuration of bodies; 
in so doing, the same grid is likely to be used. In 
particular, this may be the case in three space di- 
mensions, especially as the grid generation for this 
case typically requires a much more substantial ef- 
fort than for two space dimensions. As can be seen 
from our previous considerations, the operators T 
depend on the geometry, which does not change as 
long as the grid remains the same. These operators 
also depend on the coefficients of system (40) or (41), 
for example, on the Mach number Mo, as well as on 
the angle of attack a. Note, the angle of attack a 
formally appears neither in (40) nor in (41) since 
we always assume that the free-stream velocity is 
aligned with the positive x direction. We, however, 
take into account the angle of attack a by rotating 
the entire computational domain, see Figure 1. 

Suppose now that all the computational experi- 
ments that we are going to carry out for some cho- 
sen geometry belong to a certain range of the pa- 
rameters involved, say a m i n < a < a ma x, Mo min < 
Mo < Mo max . Then, we can pick up several points 
q>(p) within the initially prescribed range for the an- 
gle of attack and several points within the ini- 
tially prescribed range for the Mach number and 
calculate the operator T for each resulting pair 
Mq 9 ^. Note, the same is possible for the 
triplets (cv, Mo, Re)] however, the far-field effective 
Reynolds number has been noticed to influence the 
results at a much less extent than the other two pa- 
rameters. Finally, to obtain the matrix T for any 
specific pair of values (a, Mo) (within the corre- 
sponding range), we simply interpolate each entry of 
T independently with respect to the angle of attack 
and the Mach number between the known values for 

(«V), M [ 0 q) y 

We have implemented this approach numeri- 
cally for the same transonic turbulent flow around 
RAE2822 as studied above. To simplify our task 
on the preliminary stage, we used one-dimensional 
interpolation separately for a and Mo instead of us- 
ing the two-dimensional interpolation on the mesh 
M 0 (<?) J. The results shown in Table 3 corrob- 
orate usefulness of the approach based on the in- 
terpolation of coefficients of the operators T. From 



30 


V. S. RYABEN’KII and S. V. TSYNKOV 


Table 3. Interpolation of coefficients of the operator T for RAE2822 airfoil; M 0 = 0.73; Re 0 = 6.5 • 10 6 ; 
a = 2.79°; normal grid spacing near the airfoil 0.5 • 10~ 5 . 


Domain “radius” 

3 chords 

50 chords 

Grid 

600 x 104 

640 x 128 

Type of ABC’s 

p.-v. 

(59) 

Int. Mo 

Int. a 

p.-v. 

(59) 

Si 

0.8653 

0.8591 

0.8593 

0.8587 

0.8603 

0.8593 

relative error 

0.58% 

0.02% 

0.0% 

0.07% 

0% 

0% 

C d x 10 

0.1203 

0.1263 

0.1257 

0.1252 

0.1255 

0.1260 

relative error 

4.14% 

0.24% 

0.24% 

0.6% 

0% 

0% 

Cd x 10 

0.1755 

0.1816 

0.1811 

0.1805 

0.1810 

0.1815 

relative error 

3.04% 

0.05% 

0.22% 

0.55% 

0% 

0% 


this table, we see that the accuracy of the solutions 
obtained on the basis of the interpolated matrices 
is almost not worse than the accuracy that we get 
using the genuine operator T. 

Of course, the results of Table 3 are only prelimi- 
nary, and the issue of the entry-wise interpolation of 
the matrices T with respect to a and Mo requires a 
thorough further study. In particular, the approach 
based on interpolation may have certain limitations 
on the size of the computational domain, as well as 
on the actual admissible ranges of the parameters 
involved. However, the initial results are encourag- 
ing. For the repeated computations, this approach 
can drastically decrease the overall cost of the DPM- 
based ABC’s since it requires one substantial initial 
effort for calculating the matrices T on the mesh 
Mq 9 ^ , and then the boundary conditions for 

each subsequent variant of computations (within the 
prescribed range) will come for almost no extra com- 
putational cost because the cost of interpolation it- 
self is virtually negligible. 

3.4 Low Mach Number Flows 

We finally address another interesting aspect of 
implementation of the DPM-based ABC’s. It is well- 
known that many standard explicit solvers for com- 
pressible flows encounter difficulties when directly 
applied to calculating the flows with low Mach num- 
bers. The difficulties are caused by the “different 
scales” of eigenvalues u and u ± c ( u is the flow ve- 
locity and c is the speed of sound, |w,| <C c), and 
result in the severe Courant-type limitations on the 
time step. One possible cure for this problem is 
based on the so-called local preconditioning tech- 
niques. The idea of these techniques is to change the 
time-evolving system (multiplying it by some non- 
singular matrix-preconditioner) so that the gap be- 
tween the eigenvalues is narrowed but at the same 


time the steady state remains unchanged. An ap- 
proach of this type has been recently proposed by 
Turkel, Fiterman, van Leer, and Vatsa, see 4 ®’ 47, 48^ 
and has already been implemented in practice on the 
basis of the code 32, 33, 34 . 

It, however, turns out that the standard ABC’s 

no on o A 

incorporated in the code ' 5 ~’ 0,5 ’ (characteris- 
tics/extrapolation) perform poorly for the case of 
low Mach number flows. On the other hand, bound- 
ary conditions (59) in this case demonstrate the 
same good performance as they show in the case of 
transonic flows. In Table 4, we compare numerical 
results obtained using two different types of ABC’s 
for a low Mach number turbulent flow around the 
airfoil RAE2822. One can see that as in the previous 
cases, the DPM-based ABC’s allow us to maintain 
high accuracy of computations for small computa- 
tional domains. 

According to the authors of 43, 47, 43 , their pre- 
conditioning technique, actually performs better if 
the governing equations are written with respect to 
some other equivalent set of unknowns rather than 
(tq v, p, p). As concerns the DPM-based ABC’s, we 
do not reqrite equations (40) or (41). For those cases 
presented in'Table 4, we have been able to obtain 
accurate results without any changes (except in the 
input data) in the boundary conditions algorithm. 
We, however, note that in so doing the system ma- 
trices (44) become strongly non-symmetric. The ac- 
curacy of the results from Table 4 in this case is 
probably due to the fact that we solve the differ- 
ence AP by a direct method. However, for the small 
free-stream Mach numbers the use of the symmetriz- 
ers for the system matrices (for example, those pre- 
sented in work 49 ) may still be recommended. More- 
over, in our work we have constructed the three- 
dimensional DPM-based ABC’s for the true incom- 
pressible case and then implemented these boundary 
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Table 4. Low Mach number turbulent flow around RAE2822 airfoil; M 0 = 0.01; Re 0 = 6.5 • 10 6 ; a — 2.79°; 
normal grid spacing near the airfoil 0.6 • 10~ 5 . 


Domain “radius” 

2.5 chords 

20 chords 

Grid 

320 x 64 

320 x 64 

Type of ABC’s 

p.-v. 

(59) 

p.-v. 

(59) 

Ci 

0.5708 

0.5387 

0.5419 

0.5390 

relative error 

5.3% 

0.05% 

0% 

0% 

c d 

-0.0005 

0.0016 

0.0014 

0.0016 

relative error 

??? 

0% 

0% 

0% 

C D x 100 

0.5676 

0.7761 

0.7551 

0.7764 

relative error 

24.8% 

0.03% 

0% 

0% 


conditions along with the compressible solver for a 
very low Mach number (M 0 = 0.01). In this case, 
the DPM-based boundary conditions performed as 
well as they do for the higher subsonic and transonic 
Mach numbers (see * 9 ) . 

3.5 Three-Dimensional Flows 

So far, we have been mostly describing the two- 
dimensional algorithms and the two-dimensional re- 
sults. In our work 18 > 19 (see also 20 ' 21 ), we have 
constructed and implemented the nonlocal DPM- 
based ABC’s for three-dimensional steady-state vis- 
cous flows (perhaps, the most important case from 
the standpoint of current computational practice). 

As can be seen from Section 2, the ABC’s algo- 
rithm basically consists of two parts. The first part 
may be called geometrical, it comprises the construc- 
tion of the grid sets Mi n , M ex , Af in , Af ex , j, collo- 
cation grid cj, all the necessary interpolation oper- 
ations (Rrai) R W y, an d continuation 7r 7w . 

The second part is computational, it consists of the 
cross-stream transforms (e.g., discrete Fourier) and 
the solution of systems (46), (48). 

The principle differences between the two- 
dimensional case and the three-dimensional case are 
mostly concentrated in the first (geometrical) part, 
whereas the second (numerical) part changes only 
quantitatively. The computational geometry is ob- 
viously more cumbersome for the case of three space 
dimensions than it is for the case of two space dimen- 
sions. Moreover, there is a qualitative difference be- 
tween interpolation along the one-dimensional curve, 
which constitutes the artificial boundary in two di- 
mensions and which does not differ from the straight 
line when it comes to the issue of internal geom- 
etry, and interpolation along the two-dimensional 
curvilinear surface, which constitutes the artificial 
boundary in three dimensions and which has a non- 


trivial internal geometry. In practice, our compu- 
tations 18, 19 show that although the geometrical 
part of the algorithm for three space dimensions is 
more complicated, it is still universal in the sense 
that the same procedure serves a variety of artifi- 
cial boundaries with different shapes; moreover, this 
geometrical part also turns out numerically cheap. 

The major difference between the two- and three- 
dimensional cases for the computational part of the 
algorithm is that we basically add another cross- 
stream direction (in the case of a three-dimensional 
wing, for example, it may be a span-wise direction). 
Then, we separate the variables in the AP by trans- 
forming in both cross-stream and span-wise direc- 
tions and obtain a family of one-dimensional systems 
of type (46) with boundary conditions (48) that for- 
mally looks exactly the same although the matrices 
are of order five and depend on a pair of parameters 
(wavenumbers) rather than on only one parameter 
k. The numerical part of the ABC’s algorithm in 
three space dimensions also appears universal in the 
sense that it does not depend on the shape of the 
specific computational domain. The details of the 
three-dimensional DPM-based algorithm, as well as 
various computational results, can be found in our 
recent paper 19 . Here, we reproduce only one nu- 
merical example. 

We consider a steady-state flow of viscous com- 
pressible gas around the ONERA M6 wing. We use 
the NASA-developed code by Vatsa, et al. to inte- 
grate the thin-layer equations on a one-block curvi- 
linear C-0 type grid generated around the wing (the 
geometric setup for three space dimension is delin- 
eated in our work 19 ). The code 50 is similar to 
the one ^ ^ that we used for two-dimensional 

computations, it is based on the central-difference 
finite-volume discretization in space with the first- 
and third-order artificial dissipation. Pseudo-time 
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Table 5. Comparison with the standard characteristics-base boundary conditions for a turbulent flow around 
the ONERA M6 wing: M 0 = 0.84; Re 0 = 11.7 • 10 6 ; a = 3.06°. 


Domain “radius” 

3 root chords 

10 root chords 

Grid 

197 x 49 x 33 

209 x 57 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

Full lift, C L 

0.298T0.004 

0.2798 

0.2805 

0.2786 

Relative error 

6.24%±1.43% 

0.43% 

0% 

0% 

Full drag, C'd x 10 

0.168±0.008 

0.1537 

0.1542 

0.1531 

Relative error 

8.95%±5.19% 

0.39% 

0% 

0% 


iterations are used for obtaining the steady-state 
solution; the integration in time is done by the 
five-stage Runge-Kutta algorithm (with the Courant 
number calculated locally) supplemented by the 
residual smoothing. For the purpose of accelerating 
the convergence, the multigrid methodology is im- 
plemented; in our computations we used three sub- 
sequent grid levels with V cycles; the full multigrid 
methodology (FMG) could be employed as well. In 
addition, we use the preconditioning technique 5 * to 
improve the convergence to steady state. We im- 
plement the DPM-based ABC’s on the finest level 
of multigrid on the final FMG stage; the boundary 
data for coarser levels are provided by the coarsen- 
ing procedure. As mentioned above, the implemen- 
tation of the ABC’s was direct, i.e., it did not require 
the calculation of the matrices T. Moreover, even on 
the finest level we implement the DPM-based ABC’s 
only on the first and the last Runge-Kutta stages, 
which seems to make very little difference compared 
to the implementation on all five stages; the bound- 
ary data for the three intermediate stages are pro- 
vided from the DPM-based ABC’s on the first stage. 
Unlike the two-dimensional case, the standard treat- 
ment of the external boundary in three dimensions is 
based on merely the locally one- dimensional charac- 
teristics analysis and extrapolation (the point-vortex 
model is not applicable). 

We have conducted the computations for a stan- 
dard three-dimensional transonic test case for the 
ONERA M6 wing: Mo = 0.84; Re 0 = 11.7 • 

10 6 ; a = 3.06°. The solution was calculated for 
two different computational domains of the aver- 
age radii of approximately 10 and 3 root chords 
of the wing, respectively. The results summarized 
in Table 5 clearly demonstrate that for the small 
computational domains the DPM-based ABC’s gen- 
erate much more accurate solutions that the stan- 
dard (characteristics-based) boundary conditions 


do. Note, the grids for the different domains have 
different dimensions, and the smaller 197 x 49 x 33 
grid (3 root chords) is now an exact subset of the 
bigger 203 x 57 x 33 grid (10 root chords). This 
is done in order to eliminate any influence that the 
change of the grid in the near field could possibly 
exert on the solution. 

Besides the improvement of accuracy, the appli- 
cation of the DPM-based ABC’s to transonic flow 
computations on the small (3 root chords) compu- 
tational domain yielded much higher convergence 
rate of the residual (continuity equation), as well 
as much faster convergence of other quantities, in- 
cluding those deemed as sensitive, e.g., the number 
of supersonic points in the domain. In Figures 10a. 
and 10b, we show the convergence history for this su- 
percritical flow variant. One can see that the conver- 
gence for the standard boundary conditions is poor; 
therefore the corresponding force coefficients in Ta- 
ble 5 are given with the error bands indicated. 

For the 10 root chords domain, the DPM-based 
ABC’s also provide for some convergence speedup, 
although the difference between the two ABC’s tech- 
niques is less dramatic here. This is reasonable be- 
cause one could generally expect that the bigger the 
computational domain, the smaller is the influence 
that the external boundary conditions exert on the 
numerical procedure. The convergence history for 
the 10 root chords computations is shown in Fig- 
ures 11a and 11b. 

Note, from Figure 10b one can conclude that on 
the small domain the two algorithms converge to 
quite different solutions, whereas Figure lib allows 
one to assume that on the big domain the final solu- 
tions are close to one another. The data from Table 5 
corroborate these conclusions. This behavior again 
fits into the aforementioned concept that the impact 
of the ABC’s decreases as the domain size increases. 

As concerns the computational cost of the three- 
dimensional DPM-based ABC’s, we note that by ap- 
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Work 

Figure 10a. ONERA M6: Mq = 0.84, Re 0 = 11.7 - 
10 6 , a = 3.06°. Convergence history for the residual 
of the continuity equation. Average domain “radius” 
is 3 root chords of the wing; grid 197 x 49 x 33. 



Figure 11a. ONERA M6: Mq = 0.84, Re o = 11.7 • 
10 6 , a = 3.06°. Convergence history for the residual 
of the continuity equation. Average domain “radius” 
is 10 root chords of the wing; grid 209 x 57 x 33. 
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Figure 10b. ONERA M6: M 0 = 0.84, Re 0 = 11.7 • 
10°, a = 3.06°. Convergence history for the number 
of supersonic nodes in the domain. Average domain 
“radius” is 3 root chords of the wing; grid 197 X 49 x 
33. 
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Figure lib. ONERA M6: Mq = 0.84, Re o — 11.7 • 
10®, cr = 3.06°. Convergence history for the number 
of supersonic nodes in the domain. Average domain 
“radius” is 10 root chords of the wing; grid 209 x 
57 x 33. 


plying this procedure only on the first and the last 
Runge-Kutta stages and only on the finest multigrid 
level, the total number of the required calculations 
of generalized potential has been brought to a min- 
imum. In so doing, the average cost of application 
of the DPM-based ABC’s adds about 20-25% of the 
CPU time to the cost of the same procedure with 
the standard (characteristics-based) boundary con- 
ditions. This extra expense is not high (taking into 
account the improvement of accuracy); moreover, it 
can often be compensated for and even noticeably 
prevailed over by the convergence acceleration and 
the reduction of the domain size. Besides, to ex- 


plicitly decrease the computational cost associated 
with the DPM-based ABC’s the entry-wise interpo- 
lation of boundary operators (see above) and/or the 
multiresolution-based methodologies (see ^) can 
be used. We expect that the latter can also be em- 
ployed when implementing the DPM-based ABC’s 
for multi-block grids. 

4 Concluding Remarks 

The DPM-based approach provides for a geomet- 
rically universal and robust means to set the ABC’s 
for steady-state external flow computations. These 
ABC’s enable one to essentially decrease the size of 
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the computational domain (for the case of steady- 
state viscous flows) in comparison with the domains 
that can be used with standard boundary conditions. 
This implies the possibility to improve the accuracy 
of computations without increasing their cost or to 
reduce the total cost of computations without de- 
creasing the accuracy. Moreover, the total compu- 
tational cost may also be essentially reduced because 
of the convergence speedup that is accounted for by 
usage of the DPM-based ABC’s. We also emphasize 
that the DPM-based ABC’s apply to the computa- 
tional domains of irregular shape with equal ease and 
that the practical implementation of these bound- 
ary conditions is easy from the algorithmic stand- 
point. In particular, it takes a relatively little ef- 
fort to supplement some already existing code by the 
new DPM-based ABC’s. Altogether, these proper- 
ties make the DPM-based ABC’s a very attractive 
tool for external flow computations. 

The next big challenge would, of course, be ex- 
tending the DPM-based approach to the case of 
time-dependent problems. As mentioned above, a 
particular class of such problems, namely, the fluid 
flows oscillating in time, has been studied in . j n 
practice, this formulation originates, e.g., from the 
well-known problem of pitching airfoil. Analogously 
to the steady-state case, the ABC’s of ^ also con- 
nect the values of the solution at the penultimate 
and outermost rows of the grid (e.g., C-grid); how- 
ever, it is done for the entire time interval T equal 
to one period. To obtain the ABC’s in , we first 
approximate the data on T x [to, + X] (to is ar- 
bitrary) by some periodic in time vector-function in 
the sense of least squares. The approximant obvi- 
ously appears to be the Fourier series of the actual 
data on T x [to, to + X]. Assuming that the pe- 
riod T is known in advance, and that there are no 
other essential time-dependent effects in the model, 
we implement the Fourier transform in time and 
end up with the family of “steady-state” systems. 
Each member of this family has, generally speaking, 
complex coefficients but can nevertheless be treated 
analogously to how it is done above, so that the re- 
sulting boundary conditions in the frequency domain 
are obtained in the form (59). Then, implementing 
the inverse Fourier transform we obtain the ABC’s 
in time domain. When the problem is discretized in 
time, both direct and inverse transforms are discrete 
as well, and the aforementioned family of “steady- 
state” systems is finite, which gives us a way to prac- 
tically calculate the ABC’s. Clearly, the DPM-based 
ABC’s of 16 appear to be nonlocal in both space and 
time. However, the nonlocality in time is limited by 
the interval X equal to one period. 


As concerns the case of general time-dependent 
flows, for which we do not do any initial assumptions 
(like periodicity) about the behavior of the solution, 
it is, of course, more complicated and more difficult 
to handle. The main obstacle here is the nonlocality 
of the ABC’s in time. Unlike the periodic case, the 
exact ABC’s for the general time-dependent prob- 
lem would formally require storing all the preceding 
information on the artificial boundary from t = 0 to 
t = t final- As the solution develops in time, such 
boundary conditions would become more and more 
expensive from the standpoints of both memory and 
computer time, which is required for processing the 
constantly growing amount of boundary data. The 
estimated high computational cost severely limits 
the possibilities of developing the exact ABC’s for 
time-dependent problems and makes most of the 
available constructions of such boundary conditions 
practically infeasible for any long-term runs. 

Clearly, the crucial issue for any time- dependent 
ABC’s algorithm is how to effectively restrict the 
nonlocality of the boundary conditions in time. (In 
time-periodic formulation of - 16 , this restriction was 
incorporated in the formulation of the problem from 
the very beginning.) A promising approach based on 
implementation of the Laplace transform in time and 
then on usage of special recursion relations for calcu- 
lating the convolutions with the kernels of nonlocal 
operators has recently been proposed by Sofronov 
in 52 . The limitation of this technique is the re- 
quirement that the boundary should be of some reg- 
ular, e.g., linear, shape. There are several other ap- 
proaches to this problem, none of which has been 
studied thoroughly yet, although each one may in 
principle appear useful. One approach is based on 
the idea of an artificial periodic formulation in time. 
It is analogous to what we do for the cross-stream 
space coordinates in the steady-state problem. Ba- 
sically, we introduce some error, which is controlled 
by the value of the period, and in so doing obtain fi- 
nite formulation of the problem which is available for 
the numerical treatment (see Section 2). Note, the 
idea of an artificial periodic formulation in time was 
earlier proposed by Lax in 53 . Other possible ap- 
proaches could be based on some properties of solu- 
tions relevant to particular classes of equations (sys- 
tems). For example, one can make use of lacunas 
that exist in the solutions of hyperbolic equations, 
or exploit the exponential decay of coefficients of the 
Green operator as the time interval increases, which 
is relevant to parabolic equations. More details on 
these and analogous approaches can be found in the 
work 10 , in which they have been first proposed, as 
well as in the review . 
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We emphasize that the issue of ABC’s for time- 
dependent problems is important not only in CFD, 
but in many other areas of scientific computing. For 
example, the problems of computational acoustics 
and aeroacoustics (Helmholtz, wave, or linearized 
Euler equations), as well as the problems of electro- 
magnetic waves propagation (Maxwell equations), 
present significant mathematical interest and also 
attract more and more attention of the practition- 
ers. Effective algorithms for handling unbounded 
domains are required for such problems in both 
single-mode (Helmholtz-type) and wide-band (as a 
rule, time domain) formulations. For the time do- 
main algorithms, the restriction of nonlocality of the 
ABC’s in time obviously acquires particular impor- 
tance and presents a major challenge for the fu- 
ture research. Among many interesting links that 
connect this problem to other areas of computa- 
tional mathematics we would like to point out one. 
Namely, one of the principle elements of the con- 
struction of highly accurate and numerically effi- 
cient DPM-based ABC’s for time-dependent prob- 
lems would be to calculate the Calderon projections 
for the schemes that are able of clear capturing the 
lacunas relevant to the solutions of hyperbolic equa- 
tions (systems). In turn, the ideas for constructing 
such schemes seem to be closely related to another 
group of ideas and techniques associated with the 
genuinely multidimensional methods. The latter are 
intensively studied now, especially in CFD, where 
these methods present a major hope for the drastic 
increase of efficiency of the multigrid solvers (see ^ ). 
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